Modeling catalytic combustion om methane using RMG-Cat

Based somewhat on the notebook to accompany the manuscript:
"Automatic generation of microkinetic mechanisms for heterogeneous catalysis" by
C. Franklin Goldsmith, School of Engineering, Brown University, franklin_goldsmith@brown.edu, and
Richard H. West, Department of Chemical Engineering, Northeastern University, r.west@northeastern.edu

As modified to prepare data for the AIChE conference presentation: 304c Incorporation of Linear Scaling Relations into Automatic Mechanism Generation for Heterogeneous Catalysis Richard H. West, Department of Chemical Engineering, Northeastern University, Boston, MA and C. Franklin Goldsmith, Engineering, Brown University, Providence, RI Tuesday, October 31, 2017 https://aiche.confex.com/aiche/2017/meetingapp.cgi/Paper/500172 https://www.slideshare.net/richardhwest/incorporation-of-linear-scaling-relations-into-automatic-mechanism-generation-for-heterogeneous-catalysis

Then further updated to model catalytic combustion of methanol.

First, we print what git commit we were on when we ran this notebook, for both the source code (RMG-Py) and the database.

In [1]:
%%bash
cd $RMGpy
pwd
git log -n1 --pretty=oneline
cd ../RMG-database
pwd
git log -n1 --pretty=oneline
/Users/rwest/Code/Cat/RMG-Py
4f7568e85b4444cceca945b6d3116ed1105070e6 Separate the Linear Scaling Relations into new function. Apply always.
/Users/rwest/Code/Cat/RMG-database
d5bf54520fc72d16d88674014df27a8fd249c3c4 Change comments (but not values) in the recommended families list.

Model generation

We start with a base input file to generate a mechanism for CH3OH plus O2. First we print the input file we'll use to generate the model.

In [2]:
%cat base/input.py
# Data sources
database(
    thermoLibraries=['surfaceThermo', 'primaryThermoLibrary', 'thermo_DFT_CCSDTF12_BAC','DFT_QCI_thermo'],
    reactionLibraries = [('Deutschmann_Ni', False)],
    seedMechanisms = [],
    kineticsDepositories = ['training'],
    kineticsFamilies = 'default',
    kineticsEstimator = 'rate rules',
    bindingEnergies = { # default values for Ni(111)
                       'C':(-5.997, 'eV/molecule'),
                       'H':(-2.778, 'eV/molecule'),
                       'O':(-4.485, 'eV/molecule'),
                       }
)

# List of species
species(
    label='X',
    reactive=True,
    structure=adjacencyList("1 X u0"),
)

species(
    label='CH4',
    reactive=True,
    structure=SMILES("[CH4]"),
)
species(
   label='O2',
   reactive=True,
   structure=adjacencyList(
       """
1 O u1 p2 c0 {2,S}
2 O u1 p2 c0 {1,S}
"""),
)

species(
    label='N2',
    reactive=False,
    structure=SMILES("N#N"),
)

species(
    label='CO2',
    reactive=True,
    structure=SMILES("O=C=O"),
)

species(
    label='H2O',
    reactive=True,
    structure=SMILES("O"),
)

species(
    label='H2',
    reactive=True,
    structure=SMILES("[H][H]"),
)

species(
    label='CO',
    reactive=True,
    structure=SMILES("[C-]#[O+]"),
)

species(
    label='C2H6',
    reactive=True,
    structure=SMILES("CC"),
)

species(
    label='CH2O',
    reactive=True,
    structure=SMILES("C=O"),
)

species(
    label='CH3',
    reactive=True,
    structure=SMILES("[CH3]"),
)

species(
    label='C3H8',
    reactive=True,
    structure=SMILES("CCC"),
)

species(
    label='H',
    reactive=True,
    structure=SMILES("[H]"),
)

species(
    label='C2H5',
    reactive=True,
    structure=SMILES("C[CH2]"),
)

species(
    label='CH3OH',
    reactive=True,
    structure=SMILES("CO"),
)

species(
    label='HCO',
    reactive=True,
    structure=SMILES("[CH]=O"),
)

species(
    label='CH3CHO',
    reactive=True,
    structure=SMILES("CC=O"),
)

species(
    label='OH',
    reactive=True,
    structure=SMILES("[OH]"),
)

species(
    label='C2H4',
    reactive=True,
    structure=SMILES("C=C"),
)

#----------
# Reaction systems
surfaceReactor(
    temperature=(1200,'K'),
    initialPressure=(1.01325, 'bar'),
    initialGasMoleFractions={
        "CH3OH": 0.11764706,
        "O2": 0.17647059,
        "N2": 0.70588235,
    },
    initialSurfaceCoverages={
        "X": 1.0,
    },
    surfaceVolumeRatio=(1.e5, 'm^-1'),
    surfaceSiteDensity=(2.9e-9, 'mol/cm^2'),
#    terminationConversion = { "CH4":0.9,},
    terminationTime=(1.0, 's'),
)

simulator(
    atol=1e-18,
    rtol=1e-12,
)

model(
    toleranceKeepInEdge=0.0,
    toleranceMoveToCore=1e-2,
    toleranceInterruptSimulation=0.1,
    maximumEdgeSpecies=100000
)

options(
    units='si',
    saveRestartPeriod=None,
    generateOutputHTML=True,
    generatePlots=False,
    saveEdgeSpecies=True,
    saveSimulationProfiles=True,
)

Then we try running it. This will take about four minutes.

In [4]:
%%bash
python $RMGpy/rmg.py base/input.py > /dev/null
tail -n12 base/RMG.log
Saving current model edge to HTML file...
Updating RMG execution statistics...
    Execution time (DD:HH:MM:SS): 00:00:03:49
    Memory used: memory usage was unable to be logged


MODEL GENERATION COMPLETED

The final model core has 70 species and 218 reactions
The final model edge has 516 species and 977 reactions

RMG execution terminated at Tue Nov 21 23:30:27 2017

There are 70 species and 218 reactions (?)

Data processing

Next we will import some libraries and set things up to start importing and analyzing the simulation results.

In [5]:
%config InlineBackend.figure_format = 'retina'
In [6]:
%matplotlib inline
from matplotlib import pyplot as plt
import matplotlib

# The default output Type 3 (Type3) fonts can't be edited in Adobe Illustrator
# but Type 42 (TrueType) fonts can be, making it easier to move labels around
# slightly to improve layout before publication.
matplotlib.rcParams['pdf.fonttype'] = 42 

# Seaborn helps make matplotlib graphics nicer
import seaborn as sns
sns.set_style("ticks")
sns.set_context("paper", font_scale=1.5, rc={"lines.linewidth": 2.0})

import os
import re
import pandas as pd
import numpy as np
import shutil
import subprocess
import multiprocessing
In [7]:
def get_last_csv_file(job_directory):
    """
    Find the CSV file from the largest model in the provided job directory.
    
    For CSV files named `simulation_1_13.csv` you want 13 to be the highest number.
    """
    solver_directory = os.path.join(job_directory,'solver')
    csv_files = sorted([f for f in os.listdir(solver_directory) if f.endswith('.csv') ],
                       key=lambda f: int(f[:-4].split('_')[2]))
    return os.path.join(solver_directory, csv_files[-1])
    
job_directory = 'base'
get_last_csv_file(job_directory)
Out[7]:
'base/solver/simulation_1_70.csv'

We will use Pandas to import the csv file

In [8]:
last_csv_file = get_last_csv_file(job_directory)
data = pd.read_csv(last_csv_file)
data
Out[8]:
Time (s) Volume (m^3) N2 Ar He Ne X(1) CH4(2) O2(3) CO2(4) ... CH3O(51) SX(146) SX(276) SX(198) SX(250) SX(303) SX(123) SX(233) SX(87) SX(205)
0 2.659587e-24 0.098469 0.705882 0.0 0.0 0.0 0.285560 5.550766e-52 1.764706e-01 -3.000770e-120 ... 2.898438e-36 3.045506e-149 -8.214884e-112 -1.045643e-91 -1.256345e-104 -3.779405e-200 2.525411e-180 -1.666835e-229 -7.558573e-123 5.532969e-159
1 7.446844e-24 0.098469 0.705882 0.0 0.0 0.0 0.285560 8.388248e-51 1.764706e-01 5.436361e-118 ... 1.750657e-35 4.803049e-147 4.257108e-110 1.033162e-90 6.510619e-103 7.608577e-197 1.773880e-177 1.937888e-226 1.395964e-120 1.147374e-156
2 1.702136e-23 0.098469 0.705882 0.0 0.0 0.0 0.285560 9.023670e-50 1.764706e-01 2.294912e-115 ... 8.428658e-35 2.294646e-144 4.468008e-108 5.698351e-89 6.833160e-101 6.623433e-193 1.022042e-173 5.535459e-221 5.818797e-118 3.139309e-153
3 3.617039e-23 0.098469 0.705882 0.0 0.0 0.0 0.285560 8.294299e-49 1.764706e-01 7.381584e-113 ... 3.681017e-34 7.792782e-142 3.567230e-106 2.242674e-87 5.455553e-99 3.851828e-189 3.411043e-170 3.260231e-216 1.872945e-115 5.010657e-150
4 7.446844e-23 0.098469 0.705882 0.0 0.0 0.0 0.285560 7.095944e-48 1.764706e-01 2.111530e-110 ... 1.536752e-33 2.286727e-139 2.540984e-104 7.923873e-86 3.886061e-97 1.869238e-185 8.816793e-167 1.417965e-211 5.359923e-113 6.326241e-147
5 1.510646e-22 0.098469 0.705882 0.0 0.0 0.0 0.285560 5.867062e-47 1.764706e-01 5.709668e-108 ... 6.278133e-33 6.260534e-137 1.714242e-102 2.662090e-84 2.621680e-95 8.322394e-182 2.022951e-163 5.338346e-207 1.449660e-110 7.172265e-144
6 3.042568e-22 0.098469 0.705882 0.0 0.0 0.0 0.285560 4.770997e-46 1.764706e-01 1.501952e-105 ... 2.537722e-32 1.656961e-134 1.126187e-100 8.726700e-83 1.722337e-93 3.552828e-178 4.382233e-160 1.873706e-202 3.813806e-108 7.722493e-141
7 6.106412e-22 0.098469 0.705882 0.0 0.0 0.0 0.285560 3.847982e-45 1.764706e-01 3.897430e-103 ... 1.020407e-31 4.312744e-132 7.302108e-99 2.826288e-81 1.116750e-91 1.485529e-174 9.228764e-157 6.353288e-198 9.897021e-106 8.107436e-138
8 1.223410e-21 0.098469 0.705882 0.0 0.0 0.0 0.285560 3.090908e-44 1.764706e-01 1.004510e-100 ... 4.092288e-31 1.113236e-129 4.703832e-97 9.098482e-80 7.193817e-90 6.147598e-171 1.916526e-153 2.117645e-193 2.550894e-103 8.405781e-135
9 2.448948e-21 0.098469 0.705882 0.0 0.0 0.0 0.285560 2.477745e-43 1.764706e-01 2.580245e-98 ... 1.639050e-30 2.861687e-127 3.020246e-95 2.920242e-78 4.619020e-88 2.531019e-167 3.952399e-150 6.998433e-189 6.552469e-101 8.661052e-132
10 4.900023e-21 0.098469 0.705882 0.0 0.0 0.0 0.285560 1.984206e-42 1.764706e-01 6.616585e-96 ... 6.560470e-30 7.341065e-125 1.936097e-93 9.358763e-77 2.960974e-86 1.039369e-163 8.122648e-147 2.303023e-184 1.680277e-98 8.896429e-129
11 9.802174e-21 0.098469 0.705882 0.0 0.0 0.0 0.285560 1.588169e-41 1.764706e-01 1.695275e-93 ... 2.625042e-29 1.881254e-122 1.240108e-91 2.997044e-75 1.896560e-84 4.262716e-160 1.666405e-143 7.562605e-180 4.305152e-96 9.124053e-126
12 1.960648e-20 0.098469 0.705882 0.0 0.0 0.0 0.285560 1.270857e-40 1.764706e-01 4.341730e-91 ... 1.050188e-28 4.818497e-120 7.939911e-90 9.594127e-74 1.214290e-82 1.747127e-156 3.415755e-140 2.480746e-175 1.102585e-93 9.350254e-123
13 3.921508e-20 0.098469 0.705882 0.0 0.0 0.0 0.285560 1.016814e-39 1.764706e-01 1.111715e-88 ... 4.201093e-28 1.233854e-117 5.082574e-88 3.070694e-72 7.773015e-81 7.158511e-153 6.998485e-137 8.133205e-171 2.823209e-91 9.578346e-120
14 7.843229e-20 0.098469 0.705882 0.0 0.0 0.0 0.285560 8.135029e-39 1.764706e-01 2.846278e-86 ... 1.680506e-27 3.159073e-115 3.253177e-86 9.827139e-71 4.975213e-79 2.932584e-149 1.433595e-133 2.665780e-166 7.228150e-89 9.810084e-117
15 1.568667e-19 0.098469 0.705882 0.0 0.0 0.0 0.285560 6.508229e-38 1.764706e-01 7.286794e-84 ... 6.722160e-27 8.087747e-113 2.082139e-84 3.144831e-69 3.184271e-77 1.201272e-145 2.936298e-130 8.736280e-162 1.850489e-86 1.004642e-113
16 3.137355e-19 0.098469 0.705882 0.0 0.0 0.0 0.285560 5.206665e-37 1.764706e-01 1.865436e-81 ... 2.688891e-26 2.070530e-110 1.332603e-82 1.006369e-67 2.037951e-75 4.920526e-142 6.013773e-127 2.862822e-157 4.737298e-84 1.028786e-110
17 6.274732e-19 0.098469 0.705882 0.0 0.0 0.0 0.285560 4.165365e-36 1.764706e-01 4.775417e-79 ... 1.075562e-25 5.300640e-108 8.528768e-81 3.220418e-66 1.304262e-73 2.015420e-138 1.231617e-123 9.380737e-153 1.212724e-81 1.053469e-107
18 1.254949e-18 0.098469 0.705882 0.0 0.0 0.0 0.285560 3.332305e-35 1.764706e-01 1.222431e-76 ... 4.302259e-25 1.356974e-105 5.458448e-79 1.030539e-64 8.346771e-72 8.254685e-135 2.522232e-120 3.073622e-148 3.104387e-79 1.078696e-104
19 2.509899e-18 0.098469 0.705882 0.0 0.0 0.0 0.285560 2.665848e-34 1.764706e-01 3.129004e-74 ... 1.720906e-24 3.473863e-103 3.493420e-77 3.297727e-63 5.341236e-70 3.380679e-131 2.371121e-116 1.749979e-142 7.946197e-77 2.949734e-101
20 5.019801e-18 0.098469 0.705882 0.0 0.0 0.0 0.285560 2.132680e-33 1.764706e-01 8.008066e-72 ... 6.883627e-24 8.893084e-101 2.235795e-75 1.055270e-61 3.417483e-68 1.384356e-127 4.856631e-113 5.731453e-138 2.033688e-74 3.021574e-98
21 1.003960e-17 0.098469 0.705882 0.0 0.0 0.0 0.285560 1.706143e-32 1.764706e-01 2.048943e-69 ... 2.753452e-23 2.276620e-98 1.430914e-73 3.376843e-60 2.186023e-66 5.667268e-124 9.941107e-110 1.876083e-133 5.203471e-72 3.092446e-95
22 2.007921e-17 0.098469 0.705882 0.0 0.0 0.0 0.285560 1.364912e-31 1.764706e-01 5.239548e-67 ... 1.101381e-22 5.828091e-96 9.157898e-72 1.080574e-58 1.397562e-64 2.318811e-120 2.033778e-106 6.134457e-129 1.330671e-69 3.163290e-92
23 3.011881e-17 0.098469 0.705882 0.0 0.0 0.0 0.285560 4.052074e-31 1.764706e-01 3.977051e-66 ... 2.340434e-22 4.471629e-95 5.226004e-71 5.169632e-58 7.972531e-64 2.679590e-119 1.134793e-105 3.227929e-128 1.007226e-68 1.901276e-91
24 4.015842e-17 0.098469 0.705882 0.0 0.0 0.0 0.285560 8.829240e-31 1.764706e-01 1.956381e-65 ... 3.992506e-22 2.347741e-94 2.044560e-70 1.710962e-57 3.117883e-63 2.119558e-118 1.528775e-104 9.249952e-127 5.070733e-68 1.896139e-90
25 6.023763e-17 0.098469 0.705882 0.0 0.0 0.0 0.285560 2.562846e-30 1.764706e-01 2.482671e-64 ... 8.265977e-22 3.984892e-93 1.712316e-69 9.982850e-57 2.608369e-62 9.799278e-117 1.263060e-102 3.395889e-124 6.953826e-67 9.182620e-89
26 8.031684e-17 0.098469 0.705882 0.0 0.0 0.0 0.285560 5.604638e-30 1.764706e-01 1.269641e-63 ... 1.405826e-21 2.142298e-92 6.729606e-69 3.321393e-56 1.024531e-61 8.395371e-116 1.008908e-101 4.479775e-123 3.416643e-66 6.399979e-88
27 1.003961e-16 0.098469 0.705882 0.0 0.0 0.0 0.285560 1.051504e-29 1.764706e-01 4.975969e-63 ... 2.147785e-21 8.983934e-92 2.113200e-68 8.955630e-56 3.215223e-61 5.528779e-115 7.139219e-101 5.394558e-122 1.359703e-65 3.829303e-87
28 1.405545e-16 0.098469 0.705882 0.0 0.0 0.0 0.285560 2.877088e-29 1.764706e-01 1.165412e-61 ... 4.130114e-21 2.394535e-90 1.923599e-67 5.279505e-55 2.918093e-60 1.957278e-112 6.496645e-99 4.698128e-119 3.163474e-64 2.193521e-85
29 1.807129e-16 0.098469 0.705882 0.0 0.0 0.0 0.285560 6.065219e-29 1.764706e-01 7.022635e-61 ... 6.773376e-21 1.448145e-89 7.984640e-67 1.780673e-54 1.209276e-59 2.013615e-111 5.825554e-98 5.484548e-118 1.877377e-63 1.771828e-84
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
10791 9.290530e-01 0.098469 0.705882 0.0 0.0 0.0 0.239834 2.266499e-10 5.185234e-13 9.980006e-02 ... 4.428470e-12 4.895098e-15 7.794299e-12 2.649409e-14 1.578053e-10 3.945512e-17 1.781868e-09 6.599429e-13 1.659864e-09 6.645575e-11
10792 9.321807e-01 0.098469 0.705882 0.0 0.0 0.0 0.239878 2.251344e-10 5.174632e-13 9.981831e-02 ... 4.413913e-12 4.879251e-15 7.693058e-12 2.646671e-14 1.557347e-10 3.886992e-17 1.781760e-09 6.522266e-13 1.663433e-09 6.653051e-11
10793 9.353083e-01 0.098469 0.705882 0.0 0.0 0.0 0.239921 2.236381e-10 5.164090e-13 9.983648e-02 ... 4.399503e-12 4.863494e-15 7.593143e-12 2.643945e-14 1.536910e-10 3.829364e-17 1.781646e-09 6.445974e-13 1.666997e-09 6.660511e-11
10794 9.384359e-01 0.098469 0.705882 0.0 0.0 0.0 0.239965 2.221607e-10 5.153605e-13 9.985459e-02 ... 4.385240e-12 4.847828e-15 7.494537e-12 2.641231e-14 1.516738e-10 3.772613e-17 1.781528e-09 6.370544e-13 1.670557e-09 6.667955e-11
10795 9.415636e-01 0.098469 0.705882 0.0 0.0 0.0 0.240008 2.207019e-10 5.143178e-13 9.987263e-02 ... 4.371122e-12 4.832250e-15 7.397222e-12 2.638529e-14 1.496827e-10 3.716725e-17 1.781405e-09 6.295968e-13 1.674112e-09 6.675382e-11
10796 9.446912e-01 0.098469 0.705882 0.0 0.0 0.0 0.240051 2.192616e-10 5.132809e-13 9.989060e-02 ... 4.357147e-12 4.816760e-15 7.301180e-12 2.635839e-14 1.477173e-10 3.661687e-17 1.781277e-09 6.222236e-13 1.677662e-09 6.682792e-11
10797 9.478189e-01 0.098469 0.705882 0.0 0.0 0.0 0.240094 2.178395e-10 5.122495e-13 9.990850e-02 ... 4.343314e-12 4.801358e-15 7.206396e-12 2.633161e-14 1.457775e-10 3.607485e-17 1.781144e-09 6.149339e-13 1.681209e-09 6.690187e-11
10798 9.509465e-01 0.098469 0.705882 0.0 0.0 0.0 0.240137 2.164355e-10 5.112238e-13 9.992633e-02 ... 4.329621e-12 4.786042e-15 7.112852e-12 2.630494e-14 1.438628e-10 3.554107e-17 1.781007e-09 6.077269e-13 1.684750e-09 6.697565e-11
10799 9.540741e-01 0.098469 0.705882 0.0 0.0 0.0 0.240180 2.150492e-10 5.102036e-13 9.994410e-02 ... 4.316066e-12 4.770812e-15 7.020533e-12 2.627838e-14 1.419729e-10 3.501538e-17 1.780864e-09 6.006015e-13 1.688288e-09 6.704928e-11
10800 9.568890e-01 0.098469 0.705882 0.0 0.0 0.0 0.240218 2.138166e-10 5.092901e-13 9.996003e-02 ... 4.303984e-12 4.757177e-15 6.938478e-12 2.625458e-14 1.402930e-10 3.454909e-17 1.780733e-09 5.942579e-13 1.691468e-09 6.711540e-11
10801 9.597039e-01 0.098469 0.705882 0.0 0.0 0.0 0.240256 2.125980e-10 5.083810e-13 9.997591e-02 ... 4.292012e-12 4.743611e-15 6.857390e-12 2.623087e-14 1.386327e-10 3.408916e-17 1.780597e-09 5.879791e-13 1.694644e-09 6.718140e-11
10802 9.625188e-01 0.098469 0.705882 0.0 0.0 0.0 0.240294 2.113934e-10 5.074763e-13 9.999174e-02 ... 4.280148e-12 4.730112e-15 6.777257e-12 2.620725e-14 1.369918e-10 3.363551e-17 1.780458e-09 5.817646e-13 1.697817e-09 6.724727e-11
10803 9.653336e-01 0.098469 0.705882 0.0 0.0 0.0 0.240332 2.102026e-10 5.065759e-13 1.000075e-01 ... 4.268392e-12 4.716680e-15 6.698067e-12 2.618372e-14 1.353701e-10 3.318804e-17 1.780315e-09 5.756137e-13 1.700987e-09 6.731302e-11
10804 9.681485e-01 0.098469 0.705882 0.0 0.0 0.0 0.240370 2.090254e-10 5.056798e-13 1.000232e-01 ... 4.256742e-12 4.703314e-15 6.619810e-12 2.616027e-14 1.337673e-10 3.274668e-17 1.780168e-09 5.695258e-13 1.704153e-09 6.737864e-11
10805 9.709634e-01 0.098469 0.705882 0.0 0.0 0.0 0.240408 2.078618e-10 5.047880e-13 1.000389e-01 ... 4.245198e-12 4.690015e-15 6.542474e-12 2.613692e-14 1.321833e-10 3.231133e-17 1.780018e-09 5.635002e-13 1.707316e-09 6.744414e-11
10806 9.737782e-01 0.098469 0.705882 0.0 0.0 0.0 0.240445 2.067114e-10 5.039003e-13 1.000545e-01 ... 4.233759e-12 4.676781e-15 6.466049e-12 2.611365e-14 1.306178e-10 3.188191e-17 1.779864e-09 5.575365e-13 1.710475e-09 6.750952e-11
10807 9.763116e-01 0.098469 0.705882 0.0 0.0 0.0 0.240479 2.056874e-10 5.031050e-13 1.000686e-01 ... 4.223551e-12 4.664926e-15 6.398036e-12 2.609279e-14 1.292246e-10 3.150043e-17 1.779723e-09 5.522215e-13 1.713316e-09 6.756825e-11
10808 9.788450e-01 0.098469 0.705882 0.0 0.0 0.0 0.240512 2.046740e-10 5.023131e-13 1.000825e-01 ... 4.213427e-12 4.653123e-15 6.330743e-12 2.607199e-14 1.278460e-10 3.112363e-17 1.779579e-09 5.469556e-13 1.716154e-09 6.762689e-11
10809 9.811251e-01 0.098469 0.705882 0.0 0.0 0.0 0.240543 2.037708e-10 5.016032e-13 1.000951e-01 ... 4.204385e-12 4.642544e-15 6.270789e-12 2.605334e-14 1.266177e-10 3.078845e-17 1.779446e-09 5.422579e-13 1.718706e-09 6.767957e-11
10810 9.831771e-01 0.098469 0.705882 0.0 0.0 0.0 0.240570 2.029652e-10 5.009666e-13 1.001063e-01 ... 4.196304e-12 4.633059e-15 6.217320e-12 2.603660e-14 1.255222e-10 3.048995e-17 1.779325e-09 5.380634e-13 1.721001e-09 6.772692e-11
10811 9.852291e-01 0.098469 0.705882 0.0 0.0 0.0 0.240597 2.021664e-10 5.003322e-13 1.001176e-01 ... 4.188276e-12 4.623608e-15 6.164310e-12 2.601990e-14 1.244360e-10 3.019442e-17 1.779203e-09 5.339005e-13 1.723294e-09 6.777421e-11
10812 9.872812e-01 0.098469 0.705882 0.0 0.0 0.0 0.240624 2.013743e-10 4.996999e-13 1.001288e-01 ... 4.180301e-12 4.614190e-15 6.111756e-12 2.600325e-14 1.233592e-10 2.990181e-17 1.779078e-09 5.297687e-13 1.725585e-09 6.782143e-11
10813 9.893332e-01 0.098469 0.705882 0.0 0.0 0.0 0.240651 2.005888e-10 4.990697e-13 1.001400e-01 ... 4.172378e-12 4.604805e-15 6.059652e-12 2.598664e-14 1.222915e-10 2.961211e-17 1.778952e-09 5.256680e-13 1.727875e-09 6.786860e-11
10814 9.913853e-01 0.098469 0.705882 0.0 0.0 0.0 0.240677 1.998099e-10 4.984417e-13 1.001511e-01 ... 4.164507e-12 4.595453e-15 6.007996e-12 2.597008e-14 1.212330e-10 2.932528e-17 1.778824e-09 5.215981e-13 1.730163e-09 6.791569e-11
10815 9.934373e-01 0.098469 0.705882 0.0 0.0 0.0 0.240704 1.990376e-10 4.978157e-13 1.001622e-01 ... 4.156687e-12 4.586133e-15 5.956784e-12 2.595356e-14 1.201835e-10 2.904129e-17 1.778694e-09 5.175588e-13 1.732450e-09 6.796273e-11
10816 9.954894e-01 0.098469 0.705882 0.0 0.0 0.0 0.240731 1.982718e-10 4.971919e-13 1.001734e-01 ... 4.148918e-12 4.576847e-15 5.906011e-12 2.593708e-14 1.191431e-10 2.876012e-17 1.778562e-09 5.135499e-13 1.734734e-09 6.800970e-11
10817 9.975414e-01 0.098469 0.705882 0.0 0.0 0.0 0.240758 1.975125e-10 4.965701e-13 1.001844e-01 ... 4.141201e-12 4.567593e-15 5.855674e-12 2.592065e-14 1.181115e-10 2.848173e-17 1.778428e-09 5.095711e-13 1.737017e-09 6.805661e-11
10818 9.995934e-01 0.098469 0.705882 0.0 0.0 0.0 0.240784 1.967595e-10 4.959504e-13 1.001955e-01 ... 4.133533e-12 4.558371e-15 5.805770e-12 2.590426e-14 1.170887e-10 2.820610e-17 1.778293e-09 5.056222e-13 1.739299e-09 6.810346e-11
10819 1.000000e+00 0.098469 0.705882 0.0 0.0 0.0 0.240790 1.966111e-10 4.958279e-13 1.001977e-01 ... 4.132020e-12 4.556547e-15 5.795933e-12 2.590102e-14 1.168871e-10 2.815182e-17 1.778266e-09 5.048434e-13 1.739750e-09 6.811274e-11
10820 1.001645e+00 0.098469 0.705882 0.0 0.0 0.0 0.240811 1.960129e-10 4.953328e-13 1.002065e-01 ... 4.125916e-12 4.549181e-15 5.756293e-12 2.588791e-14 1.160747e-10 2.793320e-17 1.778156e-09 5.017030e-13 1.741578e-09 6.815025e-11

10821 rows × 72 columns

In [9]:
def get_pandas_data(job_directory):
    """
    Get the last CSV file from the provided job directory,
    import it into a Pandas data frame, and tidy it up a bit.
    """
    last_csv_file = get_last_csv_file(job_directory)
    data = pd.read_csv(last_csv_file)
    
    # Make the Time into an index and remove it as a column
    data.index = data['Time (s)']
    del data['Time (s)']
    # Remove inerts that RMG added automatically but we're not using
    for i in 'Ar He Ne'.split():
        del data[i]
    # Remove the Volume column
    del data['Volume (m^3)']
    # Set any amounts that are slightly negative (within the ODE solver's ATOL tolerance), equal to zero
    # to allow 'area' plots without warnings or errors.
    # Anything more negative than -1e-16 probably indicates a bug in RMG and should not be hidden here ;-)
    data[(data<0) & (data>-1e-16)] = 0
    return data
In [10]:
def rename_columns(data):
    """
    Removes the number (##) from the end of the column names, in place,
    unless doing so would make the names collide.
    Also renames a few species so the plot labels match the names in the manuscript.
    """
    import re
    old = data.columns
    new = [re.sub('\(\d+\)$','',n) for n in old]
    # don't translate names that would no longer be unique
    mapping = {k:v for k,v in zip(old,new) if new.count(v)==1}
    data.rename(columns=mapping, inplace=True)
    
    # Now change a few species that are named differently in the manuscript
    # than in the thermodynamics database used to build the model,
    # so that the plot labels match the manuscript.
    mapping = {'COX':'COvdwX', 'OCX': 'CO=X', 'C2H3X':'CH3CX', 'C2H3OX':'CH3CXO'}
    data.rename(columns=mapping, inplace=True)
In [11]:
data1 = get_pandas_data('base')
rename_columns(data1)
data1.columns
Out[11]:
Index([u'N2', u'X', u'CH4', u'O2', u'CO2', u'H2O', u'H2', u'CO', u'C2H6',
       u'CH2O', u'CH3', u'C3H8', u'H', u'C2H5', u'CH3OH', u'HCO', u'CH3CHO',
       u'OH', u'C2H4', u'CH4OX', u'OX', u'HX', u'CH3OX(39)', u'HOX', u'OCCO',
       u'CH2OX', u'CHOX', u'H2X', u'H2OX', u'CH3X', u'CH2X', u'CHX', u'CXHO',
       u'CO=X', u'COvdwX', u'CH3OX(40)', u'HO2X', u'CH4X', u'CH2O2X',
       u'CH3O2X(72)', u'CH3O2X(89)', u'CX', u'SX(68)', u'CO2X', u'SX(114)',
       u'SX(93)', u'SX(141)', u'C2H4OX', u'SX(166)', u'SX(138)', u'CH3O2X(61)',
       u'SX(169)', u'SX(98)', u'C2H5OX', u'SX(92)', u'SX(224)', u'HOCXO',
       u'CH3O', u'SX(146)', u'SX(276)', u'SX(198)', u'SX(250)', u'SX(303)',
       u'SX(123)', u'SX(233)', u'SX(87)', u'SX(205)'],
      dtype='object')
In [13]:
# Test it with some plots
data1[['CH3OH', 'O2']].plot.line()
data1[['CO', 'H2']].plot.line()
data1[['H2O']].plot.line()
Out[13]:
<matplotlib.axes._subplots.AxesSubplot at 0x10db76810>
In [15]:
species_names = data1.columns
# just the gas phase species that aren't always zero:
gas_phase = [n for n in species_names if 'X' not in n and (data1[n]>0).any()]
# all the surface species
surface_phase = [n for n in species_names if 'X' in n]
surface_phase.remove('X')
data1[gas_phase].plot.line()
data1[surface_phase].plot.line()
Out[15]:
(0, 0.05)
In [16]:
print "Significant species (those that exceed 0.001 mol at some point)"
significant = [n for n in data1.columns if(data1[n]>0.001).any()]
with sns.color_palette("hls", len(significant)):
    data1[significant].plot.area(legend='reverse')
Significant species (those that exceed 0.001 mol at some point)
In [23]:
lim = 1e-4
surface = [n for n in data1.columns if 'X' in n and n!='X' and (data1[n]>lim).any() ]
print "The {} surface species that exceed {:.1e} mol at some point".format(len(surface), lim)
total_sites = max(data1['X'])
with sns.color_palette('Set3',len(surface)):
    (data1[surface]/total_sites).plot.area(legend='reverse')
    plt.ylabel('site fraction')
    plt.xlim(0,1e-2) ## notice this is much less than 1 seconde
The 12 surface species that exceed 1.0e-04 mol at some point

Effect of binding energies

In [24]:
# First, make a series of input files in separate directories

with open(os.path.join('base', 'input.py')) as infile:
    input_file = infile.read()
    
base_directory = 'binding_energies'
def directory(carbon,oxygen):
    return os.path.join(base_directory, "c{:.3f}o{:.3f}".format(carbon,oxygen))

def make_input(binding_energies):
    """
    Make an input file for the given (carbon,oxygen) tuple (or iterable) of binding energies
    and return the name of the directory in which it is saved.
    """
    carbon, oxygen = binding_energies
    output = input_file
    out_dir = directory(carbon, oxygen)
    carbon_string = "'C':({:f}, 'eV/molecule')".format(carbon)
    output = re.sub("'C':\(.*?, 'eV/molecule'\)", carbon_string, output)
    oxygen_string = "'O':({:f}, 'eV/molecule')".format(oxygen)
    output = re.sub("'O':\(.*?, 'eV/molecule'\)", oxygen_string, output)
    os.path.exists(out_dir) or os.makedirs(out_dir)
    out_file = os.path.join(out_dir, 'input.py')
    with open(out_file,'w') as outfile:
        outfile.write(output)
    shutil.copy(os.path.join('base','run.sh'), out_dir)
    return out_dir

    
print make_input((-8,-3.5))
    
binding_energies/c-8.000o-3.500
In [16]:
def run_simulation(binding_energies):
    """
    Assuming a job file already exists, run it.  This one is local.
    Takes a tuple of binding energies, (carbon, oxygen)
    """
    carbon, oxygen = binding_energies
    job_directory = directory(carbon, oxygen)
    print job_directory
    assert os.path.exists(job_directory)
    return subprocess.check_call('./run.sh', cwd=job_directory)
# a test experiment = (-8,-3.5) make_input(experiment) run_simulation(experiment)
In [25]:
# Revised range
plt.xlim(-7.5,-2)
plt.ylim(-6.5,-1.5)
Out[25]:
(-6.5, -1.5)
In [26]:
# Revised range
carbon_range = (-7.5, -2)
oxygen_range = (-6.5, -1.5)
grid_size = 9
mesh  = np.mgrid[carbon_range[0]:carbon_range[1]:grid_size*1j, oxygen_range[0]:oxygen_range[1]:grid_size*1j]
mesh
Out[26]:
array([[[-7.5   , -7.5   , -7.5   , -7.5   , -7.5   , -7.5   , -7.5   ,
         -7.5   , -7.5   ],
        [-6.8125, -6.8125, -6.8125, -6.8125, -6.8125, -6.8125, -6.8125,
         -6.8125, -6.8125],
        [-6.125 , -6.125 , -6.125 , -6.125 , -6.125 , -6.125 , -6.125 ,
         -6.125 , -6.125 ],
        [-5.4375, -5.4375, -5.4375, -5.4375, -5.4375, -5.4375, -5.4375,
         -5.4375, -5.4375],
        [-4.75  , -4.75  , -4.75  , -4.75  , -4.75  , -4.75  , -4.75  ,
         -4.75  , -4.75  ],
        [-4.0625, -4.0625, -4.0625, -4.0625, -4.0625, -4.0625, -4.0625,
         -4.0625, -4.0625],
        [-3.375 , -3.375 , -3.375 , -3.375 , -3.375 , -3.375 , -3.375 ,
         -3.375 , -3.375 ],
        [-2.6875, -2.6875, -2.6875, -2.6875, -2.6875, -2.6875, -2.6875,
         -2.6875, -2.6875],
        [-2.    , -2.    , -2.    , -2.    , -2.    , -2.    , -2.    ,
         -2.    , -2.    ]],

       [[-6.5   , -5.875 , -5.25  , -4.625 , -4.    , -3.375 , -2.75  ,
         -2.125 , -1.5   ],
        [-6.5   , -5.875 , -5.25  , -4.625 , -4.    , -3.375 , -2.75  ,
         -2.125 , -1.5   ],
        [-6.5   , -5.875 , -5.25  , -4.625 , -4.    , -3.375 , -2.75  ,
         -2.125 , -1.5   ],
        [-6.5   , -5.875 , -5.25  , -4.625 , -4.    , -3.375 , -2.75  ,
         -2.125 , -1.5   ],
        [-6.5   , -5.875 , -5.25  , -4.625 , -4.    , -3.375 , -2.75  ,
         -2.125 , -1.5   ],
        [-6.5   , -5.875 , -5.25  , -4.625 , -4.    , -3.375 , -2.75  ,
         -2.125 , -1.5   ],
        [-6.5   , -5.875 , -5.25  , -4.625 , -4.    , -3.375 , -2.75  ,
         -2.125 , -1.5   ],
        [-6.5   , -5.875 , -5.25  , -4.625 , -4.    , -3.375 , -2.75  ,
         -2.125 , -1.5   ],
        [-6.5   , -5.875 , -5.25  , -4.625 , -4.    , -3.375 , -2.75  ,
         -2.125 , -1.5   ]]])
In [27]:
experiments = mesh.reshape((2,-1)).T
experiments
Out[27]:
array([[-7.5   , -6.5   ],
       [-7.5   , -5.875 ],
       [-7.5   , -5.25  ],
       [-7.5   , -4.625 ],
       [-7.5   , -4.    ],
       [-7.5   , -3.375 ],
       [-7.5   , -2.75  ],
       [-7.5   , -2.125 ],
       [-7.5   , -1.5   ],
       [-6.8125, -6.5   ],
       [-6.8125, -5.875 ],
       [-6.8125, -5.25  ],
       [-6.8125, -4.625 ],
       [-6.8125, -4.    ],
       [-6.8125, -3.375 ],
       [-6.8125, -2.75  ],
       [-6.8125, -2.125 ],
       [-6.8125, -1.5   ],
       [-6.125 , -6.5   ],
       [-6.125 , -5.875 ],
       [-6.125 , -5.25  ],
       [-6.125 , -4.625 ],
       [-6.125 , -4.    ],
       [-6.125 , -3.375 ],
       [-6.125 , -2.75  ],
       [-6.125 , -2.125 ],
       [-6.125 , -1.5   ],
       [-5.4375, -6.5   ],
       [-5.4375, -5.875 ],
       [-5.4375, -5.25  ],
       [-5.4375, -4.625 ],
       [-5.4375, -4.    ],
       [-5.4375, -3.375 ],
       [-5.4375, -2.75  ],
       [-5.4375, -2.125 ],
       [-5.4375, -1.5   ],
       [-4.75  , -6.5   ],
       [-4.75  , -5.875 ],
       [-4.75  , -5.25  ],
       [-4.75  , -4.625 ],
       [-4.75  , -4.    ],
       [-4.75  , -3.375 ],
       [-4.75  , -2.75  ],
       [-4.75  , -2.125 ],
       [-4.75  , -1.5   ],
       [-4.0625, -6.5   ],
       [-4.0625, -5.875 ],
       [-4.0625, -5.25  ],
       [-4.0625, -4.625 ],
       [-4.0625, -4.    ],
       [-4.0625, -3.375 ],
       [-4.0625, -2.75  ],
       [-4.0625, -2.125 ],
       [-4.0625, -1.5   ],
       [-3.375 , -6.5   ],
       [-3.375 , -5.875 ],
       [-3.375 , -5.25  ],
       [-3.375 , -4.625 ],
       [-3.375 , -4.    ],
       [-3.375 , -3.375 ],
       [-3.375 , -2.75  ],
       [-3.375 , -2.125 ],
       [-3.375 , -1.5   ],
       [-2.6875, -6.5   ],
       [-2.6875, -5.875 ],
       [-2.6875, -5.25  ],
       [-2.6875, -4.625 ],
       [-2.6875, -4.    ],
       [-2.6875, -3.375 ],
       [-2.6875, -2.75  ],
       [-2.6875, -2.125 ],
       [-2.6875, -1.5   ],
       [-2.    , -6.5   ],
       [-2.    , -5.875 ],
       [-2.    , -5.25  ],
       [-2.    , -4.625 ],
       [-2.    , -4.    ],
       [-2.    , -3.375 ],
       [-2.    , -2.75  ],
       [-2.    , -2.125 ],
       [-2.    , -1.5   ]])
In [28]:
map(make_input, experiments)
Out[28]:
['binding_energies/c-7.500o-6.500',
 'binding_energies/c-7.500o-5.875',
 'binding_energies/c-7.500o-5.250',
 'binding_energies/c-7.500o-4.625',
 'binding_energies/c-7.500o-4.000',
 'binding_energies/c-7.500o-3.375',
 'binding_energies/c-7.500o-2.750',
 'binding_energies/c-7.500o-2.125',
 'binding_energies/c-7.500o-1.500',
 'binding_energies/c-6.812o-6.500',
 'binding_energies/c-6.812o-5.875',
 'binding_energies/c-6.812o-5.250',
 'binding_energies/c-6.812o-4.625',
 'binding_energies/c-6.812o-4.000',
 'binding_energies/c-6.812o-3.375',
 'binding_energies/c-6.812o-2.750',
 'binding_energies/c-6.812o-2.125',
 'binding_energies/c-6.812o-1.500',
 'binding_energies/c-6.125o-6.500',
 'binding_energies/c-6.125o-5.875',
 'binding_energies/c-6.125o-5.250',
 'binding_energies/c-6.125o-4.625',
 'binding_energies/c-6.125o-4.000',
 'binding_energies/c-6.125o-3.375',
 'binding_energies/c-6.125o-2.750',
 'binding_energies/c-6.125o-2.125',
 'binding_energies/c-6.125o-1.500',
 'binding_energies/c-5.438o-6.500',
 'binding_energies/c-5.438o-5.875',
 'binding_energies/c-5.438o-5.250',
 'binding_energies/c-5.438o-4.625',
 'binding_energies/c-5.438o-4.000',
 'binding_energies/c-5.438o-3.375',
 'binding_energies/c-5.438o-2.750',
 'binding_energies/c-5.438o-2.125',
 'binding_energies/c-5.438o-1.500',
 'binding_energies/c-4.750o-6.500',
 'binding_energies/c-4.750o-5.875',
 'binding_energies/c-4.750o-5.250',
 'binding_energies/c-4.750o-4.625',
 'binding_energies/c-4.750o-4.000',
 'binding_energies/c-4.750o-3.375',
 'binding_energies/c-4.750o-2.750',
 'binding_energies/c-4.750o-2.125',
 'binding_energies/c-4.750o-1.500',
 'binding_energies/c-4.062o-6.500',
 'binding_energies/c-4.062o-5.875',
 'binding_energies/c-4.062o-5.250',
 'binding_energies/c-4.062o-4.625',
 'binding_energies/c-4.062o-4.000',
 'binding_energies/c-4.062o-3.375',
 'binding_energies/c-4.062o-2.750',
 'binding_energies/c-4.062o-2.125',
 'binding_energies/c-4.062o-1.500',
 'binding_energies/c-3.375o-6.500',
 'binding_energies/c-3.375o-5.875',
 'binding_energies/c-3.375o-5.250',
 'binding_energies/c-3.375o-4.625',
 'binding_energies/c-3.375o-4.000',
 'binding_energies/c-3.375o-3.375',
 'binding_energies/c-3.375o-2.750',
 'binding_energies/c-3.375o-2.125',
 'binding_energies/c-3.375o-1.500',
 'binding_energies/c-2.688o-6.500',
 'binding_energies/c-2.688o-5.875',
 'binding_energies/c-2.688o-5.250',
 'binding_energies/c-2.688o-4.625',
 'binding_energies/c-2.688o-4.000',
 'binding_energies/c-2.688o-3.375',
 'binding_energies/c-2.688o-2.750',
 'binding_energies/c-2.688o-2.125',
 'binding_energies/c-2.688o-1.500',
 'binding_energies/c-2.000o-6.500',
 'binding_energies/c-2.000o-5.875',
 'binding_energies/c-2.000o-5.250',
 'binding_energies/c-2.000o-4.625',
 'binding_energies/c-2.000o-4.000',
 'binding_energies/c-2.000o-3.375',
 'binding_energies/c-2.000o-2.750',
 'binding_energies/c-2.000o-2.125',
 'binding_energies/c-2.000o-1.500']
# Now run the simulations using a pool. ## Don't run this cell unless you have a while to wait!! ### pool = multiprocessing.Pool() result = pool.map(run_simulation, experiments)
In [21]:
base_directory = 'binding_energies'
# base_directory = 'binding_energies_local'
In [22]:
def get_data(experiment):
    carbon, oxygen = experiment
    directory(carbon,oxygen)
    data = get_pandas_data(directory(carbon,oxygen))
    rename_columns(data)
    return data
In [23]:
datas = {tuple(e): get_data(e) for e in experiments}
In [24]:
datas.keys()
Out[24]:
[(-2.0, -5.875),
 (-7.5, -3.375),
 (-6.8125, -5.875),
 (-2.0, -6.5),
 (-6.8125, -2.75),
 (-4.75, -5.875),
 (-6.125, -1.5),
 (-7.5, -6.5),
 (-2.6875, -5.875),
 (-4.75, -2.75),
 (-5.4375, -6.5),
 (-2.0, -1.5),
 (-4.75, -4.0),
 (-6.125, -4.625),
 (-4.0625, -4.625),
 (-6.125, -3.375),
 (-6.8125, -3.375),
 (-4.0625, -1.5),
 (-2.6875, -4.625),
 (-6.8125, -5.25),
 (-2.6875, -4.0),
 (-4.75, -1.5),
 (-3.375, -2.125),
 (-5.4375, -5.25),
 (-4.75, -4.625),
 (-5.4375, -2.75),
 (-2.0, -4.0),
 (-5.4375, -4.0),
 (-3.375, -2.75),
 (-6.125, -6.5),
 (-4.0625, -6.5),
 (-7.5, -5.25),
 (-2.0, -2.125),
 (-3.375, -5.875),
 (-4.75, -3.375),
 (-2.6875, -2.125),
 (-5.4375, -5.875),
 (-7.5, -5.875),
 (-3.375, -5.25),
 (-6.8125, -4.625),
 (-2.6875, -6.5),
 (-6.8125, -4.0),
 (-2.0, -4.625),
 (-6.8125, -1.5),
 (-4.0625, -2.125),
 (-3.375, -3.375),
 (-6.8125, -2.125),
 (-7.5, -2.75),
 (-2.6875, -1.5),
 (-4.0625, -5.875),
 (-4.0625, -4.0),
 (-4.0625, -5.25),
 (-4.75, -6.5),
 (-2.0, -2.75),
 (-5.4375, -1.5),
 (-3.375, -4.0),
 (-3.375, -6.5),
 (-6.125, -4.0),
 (-6.125, -2.75),
 (-4.0625, -3.375),
 (-3.375, -1.5),
 (-6.8125, -6.5),
 (-2.0, -3.375),
 (-7.5, -4.625),
 (-2.0, -5.25),
 (-2.6875, -5.25),
 (-5.4375, -3.375),
 (-4.75, -2.125),
 (-5.4375, -2.125),
 (-6.125, -5.25),
 (-2.6875, -3.375),
 (-4.0625, -2.75),
 (-6.125, -5.875),
 (-5.4375, -4.625),
 (-3.375, -4.625),
 (-7.5, -2.125),
 (-2.6875, -2.75),
 (-7.5, -1.5),
 (-7.5, -4.0),
 (-6.125, -2.125),
 (-4.75, -5.25)]
In [25]:
def get_max_co2(experiment):
    data = datas[tuple(experiment)]
    return data[['CO2']].max()
highest_co2 = max([float(get_max_co2(e)) for e in experiments])
In [26]:
ax = plt.axes()
for experiment in experiments:
    print experiment
    data = get_data(experiment)
    (data[['CO2']]/highest_co2).plot.line(ax=ax)
[-7.5 -6.5]
[-7.5   -5.875]
[-7.5  -5.25]
/Users/rwest/anaconda/envs/rmg6/lib/python2.7/site-packages/matplotlib/axes/_base.py:2917: UserWarning: Attempting to set identical left==right results
in singular transformations; automatically expanding.
left=1e-15, right=1e-15
  'left=%s, right=%s') % (left, right))
[-7.5   -4.625]
[-7.5 -4. ]
[-7.5   -3.375]
[-7.5  -2.75]
[-7.5   -2.125]
[-7.5 -1.5]
[-6.8125 -6.5   ]
[-6.8125 -5.875 ]
[-6.8125 -5.25  ]
[-6.8125 -4.625 ]
[-6.8125 -4.    ]
[-6.8125 -3.375 ]
[-6.8125 -2.75  ]
[-6.8125 -2.125 ]
[-6.8125 -1.5   ]
[-6.125 -6.5  ]
[-6.125 -5.875]
[-6.125 -5.25 ]
[-6.125 -4.625]
[-6.125 -4.   ]
[-6.125 -3.375]
[-6.125 -2.75 ]
[-6.125 -2.125]
[-6.125 -1.5  ]
[-5.4375 -6.5   ]
[-5.4375 -5.875 ]
[-5.4375 -5.25  ]
[-5.4375 -4.625 ]
[-5.4375 -4.    ]
[-5.4375 -3.375 ]
[-5.4375 -2.75  ]
[-5.4375 -2.125 ]
[-5.4375 -1.5   ]
[-4.75 -6.5 ]
[-4.75  -5.875]
[-4.75 -5.25]
[-4.75  -4.625]
[-4.75 -4.  ]
[-4.75  -3.375]
[-4.75 -2.75]
[-4.75  -2.125]
[-4.75 -1.5 ]
[-4.0625 -6.5   ]
[-4.0625 -5.875 ]
[-4.0625 -5.25  ]
[-4.0625 -4.625 ]
[-4.0625 -4.    ]
[-4.0625 -3.375 ]
[-4.0625 -2.75  ]
[-4.0625 -2.125 ]
[-4.0625 -1.5   ]
[-3.375 -6.5  ]
[-3.375 -5.875]
[-3.375 -5.25 ]
[-3.375 -4.625]
[-3.375 -4.   ]
[-3.375 -3.375]
[-3.375 -2.75 ]
[-3.375 -2.125]
[-3.375 -1.5  ]
[-2.6875 -6.5   ]
[-2.6875 -5.875 ]
[-2.6875 -5.25  ]
[-2.6875 -4.625 ]
[-2.6875 -4.    ]
[-2.6875 -3.375 ]
[-2.6875 -2.75  ]
[-2.6875 -2.125 ]
[-2.6875 -1.5   ]
[-2.  -6.5]
[-2.    -5.875]
[-2.   -5.25]
[-2.    -4.625]
[-2. -4.]
[-2.    -3.375]
[-2.   -2.75]
[-2.    -2.125]
[-2.  -1.5]
In [58]:
import seaborn as sns
plt.figure(figsize=(5, 4))
num_lines = len(experiments)
sns.set_palette(sns.color_palette("coolwarm",num_lines))

ax = plt.axes()
def make_label(experiment):
    return "{:+.1f}, {:+.1f}".format(*experiment)

for experiment in experiments:
    print experiment
    data = get_data(experiment)
    times = np.array(data.index)
    normalized = data[['CO2']].values / highest_co2
    ax.plot(times, normalized, label=make_label(experiment))
    #normalized.plot.line(ax=ax, label=make_label(experiment))
    #ax.text(times[-10], normalized[-10], make_label(experiment))

#plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.xlabel('Time (s)')
plt.ylabel('Normalized CO2 concentration')
[-7.5 -6.5]
[-7.5   -5.875]
[-7.5  -5.25]
[-7.5   -4.625]
[-7.5 -4. ]
[-7.5   -3.375]
[-7.5  -2.75]
[-7.5   -2.125]
[-7.5 -1.5]
[-6.8125 -6.5   ]
[-6.8125 -5.875 ]
[-6.8125 -5.25  ]
[-6.8125 -4.625 ]
[-6.8125 -4.    ]
[-6.8125 -3.375 ]
[-6.8125 -2.75  ]
[-6.8125 -2.125 ]
[-6.8125 -1.5   ]
[-6.125 -6.5  ]
[-6.125 -5.875]
[-6.125 -5.25 ]
[-6.125 -4.625]
[-6.125 -4.   ]
[-6.125 -3.375]
[-6.125 -2.75 ]
[-6.125 -2.125]
[-6.125 -1.5  ]
[-5.4375 -6.5   ]
[-5.4375 -5.875 ]
[-5.4375 -5.25  ]
[-5.4375 -4.625 ]
[-5.4375 -4.    ]
[-5.4375 -3.375 ]
[-5.4375 -2.75  ]
[-5.4375 -2.125 ]
[-5.4375 -1.5   ]
[-4.75 -6.5 ]
[-4.75  -5.875]
[-4.75 -5.25]
[-4.75  -4.625]
[-4.75 -4.  ]
[-4.75  -3.375]
[-4.75 -2.75]
[-4.75  -2.125]
[-4.75 -1.5 ]
[-4.0625 -6.5   ]
[-4.0625 -5.875 ]
[-4.0625 -5.25  ]
[-4.0625 -4.625 ]
[-4.0625 -4.    ]
[-4.0625 -3.375 ]
[-4.0625 -2.75  ]
[-4.0625 -2.125 ]
[-4.0625 -1.5   ]
[-3.375 -6.5  ]
[-3.375 -5.875]
[-3.375 -5.25 ]
[-3.375 -4.625]
[-3.375 -4.   ]
[-3.375 -3.375]
[-3.375 -2.75 ]
[-3.375 -2.125]
[-3.375 -1.5  ]
[-2.6875 -6.5   ]
[-2.6875 -5.875 ]
[-2.6875 -5.25  ]
[-2.6875 -4.625 ]
[-2.6875 -4.    ]
[-2.6875 -3.375 ]
[-2.6875 -2.75  ]
[-2.6875 -2.125 ]
[-2.6875 -1.5   ]
[-2.  -6.5]
[-2.    -5.875]
[-2.   -5.25]
[-2.    -4.625]
[-2. -4.]
[-2.    -3.375]
[-2.   -2.75]
[-2.    -2.125]
[-2.  -1.5]
Out[58]:
<matplotlib.text.Text at 0x132a82a90>
In [28]:
import seaborn as sns
plt.figure(figsize=(5, 4))
num_lines = len(experiments)
sns.set_palette(sns.color_palette("coolwarm",num_lines))
ax = plt.axes()

def make_label(experiment):
    return "{:+.1f}, {:+.1f}".format(*experiment)

for experiment in experiments:
    print experiment
    data = get_data(experiment)
    times = np.array(data.index)
    normalized = data[['CO2']].values[:,0] / highest_co2
    ax.loglog(times, normalized, label=make_label(experiment))
    try:
        i = (np.nonzero((np.log10(times)+np.log10(normalized)) > -10))[0][0]
    except IndexError:
        i = -1
    print i, times[i], normalized[i]
    ax.text(times[i], normalized[i], make_label(experiment), rotation=45, ha='center', va='center', fontsize=12)
    
    # plt.ylim(1e-10, 1)
    # plt.xlim(1e-10, 1)
    # plt.show()
    # ax = plt.axes()

plt.ylim(1e-10, 1)
plt.xlim(1e-10, 1)
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.xlabel('Time (s)')
plt.ylabel('Normalized CO2 concentration')
[-7.5 -6.5]
-1 1e-15 0.0
[-7.5   -5.875]
-1 4.73147159596e-28 2.18580123755e-40
[-7.5  -5.25]
/Users/rwest/anaconda/envs/rmg6/lib/python2.7/site-packages/matplotlib/ticker.py:2039: UserWarning: Data has no positive values, and therefore cannot be log-scaled.
  "Data has no positive values, and therefore cannot be "
/Users/rwest/anaconda/envs/rmg6/lib/python2.7/site-packages/ipykernel/__main__.py:17: RuntimeWarning: divide by zero encountered in log10
1966 1.01758044655e-07 0.000983583838692
[-7.5   -4.625]
1943 1.478829957e-06 6.85651114894e-05
[-7.5 -4. ]
3085 1.09689115556e-05 9.13857288297e-06
[-7.5   -3.375]
2871 1.13061859244e-05 8.84471777588e-06
[-7.5  -2.75]
2921 1.13129905e-05 8.85219364941e-06
[-7.5   -2.125]
3052 1.13238768131e-05 8.86160377795e-06
[-7.5 -1.5]
3388 1.13664805844e-05 8.81584988749e-06
[-6.8125 -6.5   ]
-1 1e-15 0.0
[-6.8125 -5.875 ]
1731 1.05516049587e-08 0.00952475127423
[-6.8125 -5.25  ]
1189 7.43821697748e-08 0.00136644755545
[-6.8125 -4.625 ]
1601 1.652092286e-06 6.09533758732e-05
[-6.8125 -4.    ]
2130 8.61496228254e-06 1.16656519119e-05
[-6.8125 -3.375 ]
2256 8.73565544361e-06 1.14804244695e-05
[-6.8125 -2.75  ]
2416 8.7296204535e-06 1.14677260749e-05
[-6.8125 -2.125 ]
2500 8.74271751203e-06 1.14660997785e-05
[-6.8125 -1.5   ]
2711 1.13435169797e-05 8.84245223298e-06
[-6.125 -6.5  ]
1544 1.30170084294e-09 0.076874366689
[-6.125 -5.875]
1029 6.70988251495e-09 0.0150325465314
[-6.125 -5.25 ]
1499 1.50911552211e-07 0.000664825251381
[-6.125 -4.625]
1741 5.36569648768e-06 1.87059241872e-05
[-6.125 -4.   ]
1923 1.66401488905e-05 6.03628452463e-06
[-6.125 -3.375]
2107 1.66883444735e-05 5.99468177056e-06
[-6.125 -2.75 ]
2297 1.66929040762e-05 5.99862133633e-06
[-6.125 -2.125]
2682 1.67202803588e-05 5.99910531778e-06
[-6.125 -1.5  ]
2917 3.90136345303e-05 2.57312470918e-06
[-5.4375 -6.5   ]
868 9.00546482952e-10 0.111834227341
[-5.4375 -5.875 ]
972 3.01375791478e-08 0.00387401252963
[-5.4375 -5.25  ]
1040 1.04596934493e-06 9.62194396337e-05
[-5.4375 -4.625 ]
1854 3.17866827029e-05 3.16823731848e-06
[-5.4375 -4.    ]
1947 5.74355853858e-05 1.74985246961e-06
[-5.4375 -3.375 ]
2064 5.75228669987e-05 1.7506173527e-06
[-5.4375 -2.75  ]
2365 5.74183921591e-05 1.74451551971e-06
[-5.4375 -2.125 ]
2636 7.9755275649e-05 1.26224864033e-06
[-5.4375 -1.5   ]
2513 0.000273782782739 3.65711879871e-07
[-4.75 -6.5 ]
1322 5.52778166353e-09 0.0182413495878
[-4.75  -5.875]
1009 2.04010427743e-07 0.000491177939462
[-4.75 -5.25]
1321 7.75234237746e-06 1.32292852504e-05
[-4.75  -4.625]
1897 0.000155032330052 6.45154417096e-07
[-4.75 -4.  ]
1970 0.000206970056451 4.86048702409e-07
[-4.75  -3.375]
2179 0.00020672550761 4.84689832291e-07
[-4.75 -2.75]
2710 0.000212123149463 4.72785081971e-07
[-4.75  -2.125]
3196 0.000564809038415 1.77243926949e-07
[-4.75 -1.5 ]
3064 0.00265988745525 3.75986552907e-08
[-4.0625 -6.5   ]
1185 4.17598581368e-08 0.00261991466909
[-4.0625 -5.875 ]
1114 1.50324048213e-06 6.70672526952e-05
[-4.0625 -5.25  ]
1535 5.62314631733e-05 1.77857485564e-06
[-4.0625 -4.625 ]
2231 0.00075484626973 1.32760507957e-07
[-4.0625 -4.    ]
2461 0.000854771869785 1.17663976991e-07
[-4.0625 -3.375 ]
3261 0.000853453359584 1.17380963714e-07
[-4.0625 -2.75  ]
2982 0.00140804081623 7.11343721338e-08
[-4.0625 -2.125 ]
3696 0.00594552153095 1.68495619816e-08
[-4.0625 -1.5   ]
-1 1.00033755396 0.0
[-3.375 -6.5  ]
1096 3.01822678401e-07 0.000351940169081
[-3.375 -5.875]
1025 1.10560914986e-05 9.13489141562e-06
[-3.375 -5.25 ]
1957 0.000390146022342 2.56729108158e-07
[-3.375 -4.625]
2442 0.00376290659147 2.66183659463e-08
[-3.375 -4.   ]
2534 0.00395302522578 2.53501935007e-08
[-3.375 -3.375]
3229 0.00427075854394 2.34536479687e-08
[-3.375 -2.75 ]
3570 0.0122048366469 8.21027206708e-09
[-3.375 -2.125]
-1 1.00252982563 0.0
[-3.375 -1.5  ]
-1 1.00252982563 0.0
[-2.6875 -6.5   ]
914 2.31215984316e-06 4.99490057523e-05
[-2.6875 -5.875 ]
911 8.43172644824e-05 1.29003655944e-06
[-2.6875 -5.25  ]
1291 0.00229807032011 4.38830015249e-08
[-2.6875 -4.625 ]
2370 0.0155812486285 6.42255712595e-09
[-2.6875 -4.    ]
2854 0.0158549206292 6.37520995673e-09
[-2.6875 -3.375 ]
3140 0.0270990301292 3.69651361492e-09
[-2.6875 -2.75  ]
-1 1.0004368958 0.0
[-2.6875 -2.125 ]
-1 1.0004368958 0.0
[-2.6875 -1.5   ]
-1 1.0004368958 0.0
[-2.  -6.5]
906 1.96221240116e-05 7.84987705412e-06
[-2.    -5.875]
1130 0.00059732756552 1.69229898963e-07
[-2.   -5.25]
1363 0.0220746110302 4.56950239262e-09
[-2.    -4.625]
1627 0.0739055974921 1.36949869326e-09
[-2. -4.]
1924 0.0842394548373 1.19190699855e-09
[-2.    -3.375]
-1 1.00151694156 0.0
[-2.   -2.75]
-1 1.00151694156 0.0
[-2.    -2.125]
-1 1.00151694156 0.0
[-2.  -1.5]
-1 1.00151694156 0.0
Out[28]:
<matplotlib.text.Text at 0x11681da10>
In [29]:
ax = plt.axes()
for experiment in experiments:
    print experiment
    data = get_data(experiment)
    normalized = data[['CO2']] / highest_co2
    linearized = -np.log( 1 - normalized)
    linearized.plot.line(ax=ax)
    
[-7.5 -6.5]
[-7.5   -5.875]
[-7.5  -5.25]
[-7.5   -4.625]
[-7.5 -4. ]
[-7.5   -3.375]
[-7.5  -2.75]
[-7.5   -2.125]
[-7.5 -1.5]
[-6.8125 -6.5   ]
[-6.8125 -5.875 ]
[-6.8125 -5.25  ]
[-6.8125 -4.625 ]
[-6.8125 -4.    ]
[-6.8125 -3.375 ]
[-6.8125 -2.75  ]
[-6.8125 -2.125 ]
[-6.8125 -1.5   ]
[-6.125 -6.5  ]
[-6.125 -5.875]
[-6.125 -5.25 ]
[-6.125 -4.625]
[-6.125 -4.   ]
[-6.125 -3.375]
[-6.125 -2.75 ]
[-6.125 -2.125]
[-6.125 -1.5  ]
[-5.4375 -6.5   ]
[-5.4375 -5.875 ]
[-5.4375 -5.25  ]
/Users/rwest/anaconda/envs/rmg6/lib/python2.7/site-packages/ipykernel/__main__.py:6: RuntimeWarning: divide by zero encountered in log
[-5.4375 -4.625 ]
[-5.4375 -4.    ]
[-5.4375 -3.375 ]
[-5.4375 -2.75  ]
[-5.4375 -2.125 ]
[-5.4375 -1.5   ]
[-4.75 -6.5 ]
[-4.75  -5.875]
[-4.75 -5.25]
[-4.75  -4.625]
[-4.75 -4.  ]
[-4.75  -3.375]
[-4.75 -2.75]
[-4.75  -2.125]
[-4.75 -1.5 ]
[-4.0625 -6.5   ]
[-4.0625 -5.875 ]
[-4.0625 -5.25  ]
[-4.0625 -4.625 ]
[-4.0625 -4.    ]
[-4.0625 -3.375 ]
[-4.0625 -2.75  ]
[-4.0625 -2.125 ]
[-4.0625 -1.5   ]
[-3.375 -6.5  ]
[-3.375 -5.875]
[-3.375 -5.25 ]
[-3.375 -4.625]
[-3.375 -4.   ]
[-3.375 -3.375]
[-3.375 -2.75 ]
[-3.375 -2.125]
[-3.375 -1.5  ]
[-2.6875 -6.5   ]
[-2.6875 -5.875 ]
[-2.6875 -5.25  ]
[-2.6875 -4.625 ]
[-2.6875 -4.    ]
[-2.6875 -3.375 ]
[-2.6875 -2.75  ]
[-2.6875 -2.125 ]
[-2.6875 -1.5   ]
[-2.  -6.5]
[-2.    -5.875]
[-2.   -5.25]
[-2.    -4.625]
[-2. -4.]
[-2.    -3.375]
[-2.   -2.75]
[-2.    -2.125]
[-2.  -1.5]
In [30]:
sns.set_palette('Set1')
x_data = np.array(linearized.index)
y_data = linearized.values[:,0] 
plt.plot(x_data, y_data)
plt.show()

import scipy.stats
slope, intercept, rvalue, pvalue, stderr = scipy.stats.linregress(x_data, y_data)
print(slope, intercept, rvalue, pvalue, stderr)

print("Slope: {}".format(slope))
print("Intercept: {}".format(intercept))
print("Coefficient of determination (r squared): {}".format(rvalue*rvalue))
print("p-value (probability that the slope is zero): {}".format(pvalue))
print("Standard error in slope: {}".format(stderr))
plt.plot(x_data, y_data, '.')
plt.plot(x_data, x_data*slope+intercept)
plt.show()
(0.0, 0.0, 0.0, 1.0, 0.0)
Slope: 0.0
Intercept: 0.0
Coefficient of determination (r squared): 0.0
p-value (probability that the slope is zero): 1.0
Standard error in slope: 0.0
In [56]:
sns.set_palette('Set1')

def my_function(time, rate):
    "Thing we want to fit."
    return 1. - np.exp(-1*time*rate)
import scipy.optimize


def fit_rate(data):
    normalized = data[['CO2']] / highest_co2
    x_data = np.array(normalized.index)
    y_data = normalized.values[:,0] 

    popt, pcov = scipy.optimize.curve_fit(my_function, x_data, y_data)
    optimal_parameters = popt
    parameter_errors = np.sqrt(np.diag(pcov))
    print("Rate: {} +/- {} (1 st. dev.)".format(optimal_parameters[0],parameter_errors[0]))

    plt.plot(x_data, y_data, 'o')
    plt.plot(x_data, my_function(x_data, *optimal_parameters))
    plt.xlabel('Time (s)')
    plt.ylabel('Normalized CO2 concentration')
    
    plt.show()
    fitted_rate = optimal_parameters[0]
    return fitted_rate

rates = []
for experiment in experiments:
    print experiment
    data = get_data(experiment)
    rate = fit_rate(data)
    print rate
    rates.append(rate)

rates
    
[-7.5 -6.5]
Rate: 1.0 +/- inf (1 st. dev.)
1.0
[-7.5   -5.875]
Rate: 1.0 +/- inf (1 st. dev.)
1.0
[-7.5  -5.25]
Rate: 0.0787184234775 +/- 0.0018577897442 (1 st. dev.)
0.0787184234775
[-7.5   -4.625]
Rate: 0.0456170493524 +/- 0.00128164518463 (1 st. dev.)
0.0456170493524
[-7.5 -4. ]
Rate: 0.0499007736262 +/- 0.000610642019891 (1 st. dev.)
0.0499007736262
[-7.5   -3.375]
Rate: 0.0562935536737 +/- 0.000782533379151 (1 st. dev.)
0.0562935536737
[-7.5  -2.75]
Rate: 0.0859355700246 +/- 0.00103204695101 (1 st. dev.)
0.0859355700246
[-7.5   -2.125]
Rate: 0.0587907477575 +/- 0.000810662618911 (1 st. dev.)
0.0587907477575
[-7.5 -1.5]
Rate: 0.0456022022073 +/- 0.000434326360069 (1 st. dev.)
0.0456022022073
[-6.8125 -6.5   ]
Rate: 1.0 +/- inf (1 st. dev.)
1.0
[-6.8125 -5.875 ]
Rate: 3.86762016195 +/- 0.119274774976 (1 st. dev.)
3.86762016195
[-6.8125 -5.25  ]
Rate: 0.899774842765 +/- 0.044553038577 (1 st. dev.)
0.899774842765
[-6.8125 -4.625 ]
Rate: 1.99513984368 +/- 0.0392468255952 (1 st. dev.)
1.99513984368
[-6.8125 -4.    ]
Rate: 1.67184268389 +/- 0.0182998527649 (1 st. dev.)
1.67184268389
[-6.8125 -3.375 ]
Rate: 1.86531739755 +/- 0.0171142125846 (1 st. dev.)
1.86531739755
[-6.8125 -2.75  ]
Rate: 1.81464573046 +/- 0.0164545397689 (1 st. dev.)
1.81464573046
[-6.8125 -2.125 ]
Rate: 1.06983390676 +/- 0.00929820743948 (1 st. dev.)
1.06983390676
[-6.8125 -1.5   ]
Rate: 0.448943216684 +/- 0.00309400207326 (1 st. dev.)
0.448943216684
[-6.125 -6.5  ]
Rate: 0.618895495476 +/- 0.0449012035554 (1 st. dev.)
0.618895495476
[-6.125 -5.875]
Rate: 402021.82749 +/- 10906.2569989 (1 st. dev.)
402021.82749
[-6.125 -5.25 ]
Rate: 926.24993947 +/- 9.6661076419 (1 st. dev.)
926.24993947
[-6.125 -4.625]
Rate: 12.4885107598 +/- 0.0143531291963 (1 st. dev.)
12.4885107598
[-6.125 -4.   ]
Rate: 6.31592927129 +/- 0.0141127454062 (1 st. dev.)
6.31592927129
[-6.125 -3.375]
Rate: 5.97743867029 +/- 0.0103476947341 (1 st. dev.)
5.97743867029
[-6.125 -2.75 ]
Rate: 3.8194720879 +/- 0.013669025976 (1 st. dev.)
3.8194720879
[-6.125 -2.125]
Rate: 1.10273493857 +/- 0.00653873928067 (1 st. dev.)
1.10273493857
[-6.125 -1.5  ]
Rate: 0.194323599721 +/- 0.000516564820287 (1 st. dev.)
0.194323599721
[-5.4375 -6.5   ]
Rate: 11481231.2801 +/- 356498.395778 (1 st. dev.)
11481231.2801
[-5.4375 -5.875 ]
Rate: 8833.09869363 +/- 176.114452835 (1 st. dev.)
8833.09869363
[-5.4375 -5.25  ]
Rate: 57.1496884977 +/- 0.349095563166 (1 st. dev.)
57.1496884977
[-5.4375 -4.625 ]
Rate: 0.545452439186 +/- 0.000534438465412 (1 st. dev.)
0.545452439186
[-5.4375 -4.    ]
Rate: 0.395643853647 +/- 0.000456278881411 (1 st. dev.)
0.395643853647
[-5.4375 -3.375 ]
Rate: 0.366457840896 +/- 0.000412266459795 (1 st. dev.)
0.366457840896
[-5.4375 -2.75  ]
Rate: 0.300279482204 +/- 0.000332330375822 (1 st. dev.)
0.300279482204
[-5.4375 -2.125 ]
Rate: 0.0696365306141 +/- 0.000116743682657 (1 st. dev.)
0.0696365306141
[-5.4375 -1.5   ]
Rate: 0.00342013057508 +/- 1.63331294899e-06 (1 st. dev.)
0.00342013057508
[-4.75 -6.5 ]
Rate: 280336.305516 +/- 6798.76012658 (1 st. dev.)
280336.305516
[-4.75  -5.875]
Rate: 168.592533981 +/- 2.93151744498 (1 st. dev.)
168.592533981
[-4.75 -5.25]
Rate: 1.21785477115 +/- 0.00361577504812 (1 st. dev.)
1.21785477115
[-4.75  -4.625]
Rate: 0.0241012616654 +/- 2.12295623244e-05 (1 st. dev.)
0.0241012616654
[-4.75 -4.  ]
Rate: 0.0179200336311 +/- 4.96760018098e-06 (1 st. dev.)
0.0179200336311
[-4.75  -3.375]
Rate: 0.0179821021011 +/- 5.23913914904e-06 (1 st. dev.)
0.0179821021011
[-4.75 -2.75]
Rate: 0.0157115510532 +/- 1.18548308323e-05 (1 st. dev.)
0.0157115510532
[-4.75  -2.125]
Rate: 0.00133176500066 +/- 6.40859259773e-07 (1 st. dev.)
0.00133176500066
[-4.75 -1.5 ]
Rate: 2.750220367e-05 +/- 7.61733264403e-09 (1 st. dev.)
2.750220367e-05
[-4.0625 -6.5   ]
Rate: 3804.85264633 +/- 90.3786484269 (1 st. dev.)
3804.85264633
[-4.0625 -5.875 ]
Rate: 3.77052411673 +/- 0.0525730147841 (1 st. dev.)
3.77052411673
[-4.0625 -5.25  ]
Rate: 0.0343355785466 +/- 2.75173621853e-06 (1 st. dev.)
0.0343355785466
[-4.0625 -4.625 ]
Rate: 0.000931372424518 +/- 1.38749969894e-06 (1 st. dev.)
0.000931372424518
[-4.0625 -4.    ]
Rate: 0.000825223852457 +/- 1.12630402581e-06 (1 st. dev.)
0.000825223852457
[-4.0625 -3.375 ]
Rate: 0.000840567730918 +/- 1.10960289601e-06 (1 st. dev.)
0.000840567730918
[-4.0625 -2.75  ]
Rate: 0.000275372913018 +/- 3.52601251719e-07 (1 st. dev.)
0.000275372913018
[-4.0625 -2.125 ]
Rate: 1.32034761598e-05 +/- 1.37934110155e-08 (1 st. dev.)
1.32034761598e-05
[-4.0625 -1.5   ]
Rate: -2.33333538808e-10 +/- inf (1 st. dev.)
-2.33333538808e-10
[-3.375 -6.5  ]
Rate: 65.8048656546 +/- 1.1597241942 (1 st. dev.)
65.8048656546
[-3.375 -5.875]
Rate: 0.549252088504 +/- 0.00125277228127 (1 st. dev.)
0.549252088504
[-3.375 -5.25 ]
Rate: 0.000632739548047 +/- 6.06729168612e-08 (1 st. dev.)
0.000632739548047
[-3.375 -4.625]
Rate: 3.89379459154e-05 +/- 8.2323191566e-08 (1 st. dev.)
3.89379459154e-05
[-3.375 -4.   ]
Rate: 4.05486514707e-05 +/- 6.21647127033e-08 (1 st. dev.)
4.05486514707e-05
[-3.375 -3.375]
Rate: 3.22468823844e-05 +/- 5.54632148646e-08 (1 st. dev.)
3.22468823844e-05
[-3.375 -2.75 ]
Rate: 2.21988549271e-06 +/- 3.35678897056e-09 (1 st. dev.)
2.21988549271e-06
[-3.375 -2.125]
Rate: -1.71270658814e-10 +/- inf (1 st. dev.)
-1.71270658814e-10
[-3.375 -1.5  ]
Rate: -1.71270658814e-10 +/- inf (1 st. dev.)
-1.71270658814e-10
[-2.6875 -6.5   ]
Rate: 1.82785646918 +/- 0.0253622081972 (1 st. dev.)
1.82785646918
[-2.6875 -5.875 ]
Rate: 0.01511855174 +/- 1.08241703007e-06 (1 st. dev.)
0.01511855174
[-2.6875 -5.25  ]
Rate: 1.24047528663e-05 +/- 3.18453952441e-09 (1 st. dev.)
1.24047528663e-05
[-2.6875 -4.625 ]
Rate: 1.46131502279e-06 +/- 4.03354162912e-09 (1 st. dev.)
1.46131502279e-06
[-2.6875 -4.    ]
Rate: 1.50602647282e-06 +/- 3.73918122457e-09 (1 st. dev.)
1.50602647282e-06
[-2.6875 -3.375 ]
Rate: 6.78029839621e-07 +/- 1.98625861994e-09 (1 st. dev.)
6.78029839621e-07
[-2.6875 -2.75  ]
Rate: -1.91856987865e-10 +/- inf (1 st. dev.)
-1.91856987865e-10
[-2.6875 -2.125 ]
Rate: -1.91856987865e-10 +/- inf (1 st. dev.)
-1.91856987865e-10
[-2.6875 -1.5   ]
Rate: -1.91856987865e-10 +/- inf (1 st. dev.)
-1.91856987865e-10
[-2.  -6.5]
Rate: 0.322454743552 +/- 0.000513674257576 (1 st. dev.)
0.322454743552
[-2.    -5.875]
Rate: 0.000283302553654 +/- 3.60162770153e-10 (1 st. dev.)
0.000283302553654
[-2.   -5.25]
Rate: 2.5251129205e-07 +/- 1.82687038688e-10 (1 st. dev.)
2.5251129205e-07
[-2.    -4.625]
Rate: 5.27831360154e-08 +/- 1.83755836353e-10 (1 st. dev.)
5.27831360154e-08
[-2. -4.]
Rate: 3.6652205605e-08 +/- 1.06477340825e-10 (1 st. dev.)
3.6652205605e-08
[-2.    -3.375]
Rate: -2.28800906778e-10 +/- inf (1 st. dev.)
-2.28800906778e-10
[-2.   -2.75]
Rate: -2.28800906778e-10 +/- inf (1 st. dev.)
-2.28800906778e-10
[-2.    -2.125]
Rate: -2.28800906778e-10 +/- inf (1 st. dev.)
-2.28800906778e-10
[-2.  -1.5]
Rate: -2.28800906778e-10 +/- inf (1 st. dev.)
-2.28800906778e-10
Out[56]:
[1.0,
 1.0,
 0.078718423477541127,
 0.045617049352364519,
 0.049900773626224547,
 0.05629355367373081,
 0.085935570024571611,
 0.058790747757543366,
 0.04560220220731221,
 1.0,
 3.8676201619457853,
 0.89977484276486119,
 1.9951398436782986,
 1.6718426838861284,
 1.8653173975476403,
 1.8146457304581183,
 1.0698339067596483,
 0.44894321668425724,
 0.61889549547613942,
 402021.82748993218,
 926.24993946968596,
 12.488510759781825,
 6.3159292712937649,
 5.9774386702895423,
 3.8194720879032857,
 1.1027349385654395,
 0.19432359972074476,
 11481231.280071581,
 8833.0986936290392,
 57.149688497651304,
 0.54545243918564312,
 0.39564385364665822,
 0.36645784089578448,
 0.30027948220406198,
 0.069636530614081732,
 0.0034201305750834517,
 280336.30551596911,
 168.59253398056438,
 1.2178547711537022,
 0.02410126166543276,
 0.01792003363105989,
 0.017982102101145765,
 0.015711551053246273,
 0.001331765000661682,
 2.7502203670048636e-05,
 3804.8526463329972,
 3.770524116734403,
 0.03433557854664341,
 0.00093137242451807485,
 0.0008252238524569125,
 0.00084056773091780823,
 0.00027537291301832915,
 1.3203476159808562e-05,
 -2.3333353880815962e-10,
 65.804865654628898,
 0.54925208850394092,
 0.00063273954804698898,
 3.8937945915391896e-05,
 4.0548651470687108e-05,
 3.2246882384376838e-05,
 2.2198854927098678e-06,
 -1.7127065881364985e-10,
 -1.7127065881364985e-10,
 1.827856469182163,
 0.015118551740000472,
 1.2404752866319034e-05,
 1.4613150227874964e-06,
 1.5060264728215717e-06,
 6.7802983962091818e-07,
 -1.9185698786507375e-10,
 -1.9185698786507375e-10,
 -1.9185698786507375e-10,
 0.32245474355219145,
 0.00028330255365353977,
 2.5251129205000391e-07,
 5.2783136015371981e-08,
 3.6652205604959182e-08,
 -2.2880090677751836e-10,
 -2.2880090677751836e-10,
 -2.2880090677751836e-10,
 -2.2880090677751836e-10]
In [32]:
rates = np.array(rates)
fixed_rates = rates * (rates>0) + (1e-9 * (rates<0))
log_rates = np.log(fixed_rates)
log_rates
Out[32]:
array([  0.        ,   0.        ,  -2.54187805,  -3.08747374,
        -2.99771877,  -2.87717525,  -2.45415745,  -2.83377079,
        -3.08779927,   0.        ,   1.35263937,  -0.10561072,
         0.69071414,   0.51392642,   0.62343122,   0.59589026,
         0.06750341,  -0.80085887,  -0.47981885,  12.90426166,
         6.83114411,   2.52480908,   1.8430749 ,   1.78799216,
         1.34011222,   0.0977934 ,  -1.63823047,  16.2562242 ,
         9.08626116,   4.04567394,  -0.60613967,  -0.92724083,
        -1.0038718 ,  -1.20304163,  -2.66446598,  -5.67807655,
        12.54374525,   5.12748476,   0.19709093,  -3.72549109,
        -4.02183599,  -4.01837834,  -4.1533591 ,  -6.62125015,
       -10.50124442,   8.24403254,   1.32721401,  -3.37157319,
        -6.97885133,  -7.09985587,  -7.08143302,  -8.19738433,
       -11.23503042, -20.72326584,   4.18669378,  -0.59919777,
        -7.36545168, -10.15354131, -10.11300803, -10.34208919,
       -13.01805494, -20.72326584, -20.72326584,   0.60314395,
        -4.1918327 , -11.29743086, -13.43617383, -13.40603585,
       -14.20407454, -20.72326584, -20.72326584, -20.72326584,
        -1.13179248,  -8.16899514, -15.19180987, -16.75707409,
       -17.12179223, -20.72326584, -20.72326584, -20.72326584, -20.72326584])
In [33]:
experiments
Out[33]:
array([[-7.5   , -6.5   ],
       [-7.5   , -5.875 ],
       [-7.5   , -5.25  ],
       [-7.5   , -4.625 ],
       [-7.5   , -4.    ],
       [-7.5   , -3.375 ],
       [-7.5   , -2.75  ],
       [-7.5   , -2.125 ],
       [-7.5   , -1.5   ],
       [-6.8125, -6.5   ],
       [-6.8125, -5.875 ],
       [-6.8125, -5.25  ],
       [-6.8125, -4.625 ],
       [-6.8125, -4.    ],
       [-6.8125, -3.375 ],
       [-6.8125, -2.75  ],
       [-6.8125, -2.125 ],
       [-6.8125, -1.5   ],
       [-6.125 , -6.5   ],
       [-6.125 , -5.875 ],
       [-6.125 , -5.25  ],
       [-6.125 , -4.625 ],
       [-6.125 , -4.    ],
       [-6.125 , -3.375 ],
       [-6.125 , -2.75  ],
       [-6.125 , -2.125 ],
       [-6.125 , -1.5   ],
       [-5.4375, -6.5   ],
       [-5.4375, -5.875 ],
       [-5.4375, -5.25  ],
       [-5.4375, -4.625 ],
       [-5.4375, -4.    ],
       [-5.4375, -3.375 ],
       [-5.4375, -2.75  ],
       [-5.4375, -2.125 ],
       [-5.4375, -1.5   ],
       [-4.75  , -6.5   ],
       [-4.75  , -5.875 ],
       [-4.75  , -5.25  ],
       [-4.75  , -4.625 ],
       [-4.75  , -4.    ],
       [-4.75  , -3.375 ],
       [-4.75  , -2.75  ],
       [-4.75  , -2.125 ],
       [-4.75  , -1.5   ],
       [-4.0625, -6.5   ],
       [-4.0625, -5.875 ],
       [-4.0625, -5.25  ],
       [-4.0625, -4.625 ],
       [-4.0625, -4.    ],
       [-4.0625, -3.375 ],
       [-4.0625, -2.75  ],
       [-4.0625, -2.125 ],
       [-4.0625, -1.5   ],
       [-3.375 , -6.5   ],
       [-3.375 , -5.875 ],
       [-3.375 , -5.25  ],
       [-3.375 , -4.625 ],
       [-3.375 , -4.    ],
       [-3.375 , -3.375 ],
       [-3.375 , -2.75  ],
       [-3.375 , -2.125 ],
       [-3.375 , -1.5   ],
       [-2.6875, -6.5   ],
       [-2.6875, -5.875 ],
       [-2.6875, -5.25  ],
       [-2.6875, -4.625 ],
       [-2.6875, -4.    ],
       [-2.6875, -3.375 ],
       [-2.6875, -2.75  ],
       [-2.6875, -2.125 ],
       [-2.6875, -1.5   ],
       [-2.    , -6.5   ],
       [-2.    , -5.875 ],
       [-2.    , -5.25  ],
       [-2.    , -4.625 ],
       [-2.    , -4.    ],
       [-2.    , -3.375 ],
       [-2.    , -2.75  ],
       [-2.    , -2.125 ],
       [-2.    , -1.5   ]])
In [51]:
rate_grid = np.reshape(log_rates, (grid_size,grid_size))
In [52]:
ax = sns.heatmap(rate_grid)
In [36]:
extent = carbon_range + oxygen_range
extent
Out[36]:
(-7.5, -2, -6.5, -1.5)
In [37]:
plt.imshow(rate_grid, interpolation='spline16', origin='lower', extent=extent, aspect='equal')
plt.plot(-5.997, -4.485, 'ok')
plt.text(-5.997, -4.485, 'Ni(111)')
Out[37]:
<matplotlib.text.Text at 0x1170af050>
In [38]:
# (1)	Medford, A. J.; Lausche, A. C.; Abild-Pedersen, F.; Temel, B.; Schjødt, N. C.; Nørskov, J. K.; Studt, F. Activity and Selectivity Trends in Synthesis Gas Conversion to Higher Alcohols. Topics in Catalysis 2014, 57 (1-4), 135–142 DOI: 10.1007/s11244-013-0169-0.

medford_energies = { # Carbon, then Oxygen
'Ru': ( 0.010349288486416697, -2.8153856448231256),
'Rh': ( 0.16558861578266493, -2.546620868091181),
'Ni': ( 0.3001293661060802, -2.5881741535441853),
'Ir': ( 0.36222509702457995, -2.826185484230718),
'Pd': ( 0.28460543337645516, -1.207119596734621),
'Pt': ( 0.8796895213454077, -1.445820136503547),
'Cu': ( 2.323415265200518, -1.7218249542757729),
'Ag': ( 3.855109961190168, -0.8341504215550701),
'Au': ( 3.5601552393272975, -0.10963108355266138),
}
In [39]:
# Shift medford's energies so that Ni matches Wayne Blaylock's Ni
blaylock_ni = np.array([-5.997, -4.485])
old_ni = np.array(medford_energies['Ni'])
shifted_energies = {metal: tuple(blaylock_ni + np.array(E)-old_ni) for metal,E in medford_energies.items()}
shifted_energies
Out[39]:
{'Ag': (-2.4420194049159121, -2.7309762680108851),
 'Au': (-2.7369741267787826, -2.006456930008476),
 'Cu': (-3.9737141009055619, -3.6186508007315878),
 'Ir': (-5.9349042690814997, -4.723011330686532),
 'Ni': (-5.9969999999999999, -4.4849999999999994),
 'Pd': (-6.0125239327296249, -3.1039454431904363),
 'Pt': (-5.417439844760672, -3.3426459829593616),
 'Rh': (-6.1315407503234152, -4.4434467145469956),
 'Ru': (-6.2867800776196638, -4.712211491278941)}
In [40]:
plt.imshow(rate_grid, interpolation='spline16', origin='lower', extent=extent, aspect='equal')
for metal, coords in shifted_energies.iteritems():
    plt.plot(coords[0], coords[1], 'ok')
    plt.text(coords[0], coords[1], metal)
plt.xlim(carbon_range)
plt.ylim(oxygen_range)
Out[40]:
(-6.5, -1.5)
In [41]:
# For close packed surfaces from
# Abild-Pedersen, F.; Greeley, J.; Studt, F.; Rossmeisl, J.; Munter, T. R.; Moses, P. G.; Skúlason, E.; Bligaard, T.; Norskov, J. K. Scaling Properties of Adsorption Energies for Hydrogen-Containing Molecules on Transition-Metal Surfaces. Phys. Rev. Lett. 2007, 99 (1), 016105 DOI: 10.1103/PhysRevLett.99.016105.
abildpedersen_energies = { # Carbon, then Oxygen
'Ru': ( -6.397727272727272, -5.104763568600047),
'Rh': ( -6.5681818181818175, -4.609771721406942),
'Ni': ( -6.045454545454545, -4.711681807593758),
'Ir': ( -6.613636363636363, -5.94916142557652),
'Pd': ( -6, -3.517877940833916),
'Pt': ( -6.363636363636363, -3.481481481481482),
'Cu': ( -4.159090909090907, -3.85272536687631),
'Ag': ( -2.9545454545454533, -2.9282552993244817),
'Au': ( -3.7499999999999973, -2.302236198462614),
}
In [141]:
abildpedersen_energies['Pt'][0] - abildpedersen_energies['Ni'][0]
Out[141]:
-0.31818181818181834
In [114]:
c_step = mesh[0,1,0]-mesh[0,0,0]
o_step = mesh[1,0,1]-mesh[1,0,0]
carbon_range2 = (carbon_range[0]-c_step/2, carbon_range[1]+c_step/2)
oxygen_range2 = (oxygen_range[0]-c_step/2, oxygen_range[1]+c_step/2)
extent2 = carbon_range2 + oxygen_range2
In [139]:
plt.imshow(rate_grid.T, interpolation='spline16', origin='lower', extent=extent2, aspect='equal')
for metal, coords in abildpedersen_energies.iteritems():
    color = {'Ag':'w','Au':'w','Cu':'w'}.get(metal,'k')
    plt.plot(coords[0], coords[1], 'o'+color)
    plt.text(coords[0], coords[1], metal, color=color)
plt.xlim(carbon_range)
plt.ylim(oxygen_range)
plt.xlabel('$\Delta E^C$ (eV)')
plt.ylabel('$\Delta E^O$ (eV)')
Out[139]:
<matplotlib.text.Text at 0x133dfa990>
In [43]:
new_rates = []
for experiment in experiments:
    print experiment
    data = get_data(experiment)
    times = np.array(data.index)
    normalized = data[['CO2']].values[:,0] / highest_co2
    try:
        i = (np.nonzero((np.log10(times)+np.log10(normalized)) > -10))[0][0]
        time = times[i]
    except IndexError:
        time = 1.0
    new_rates.append(1./time)
new_log_rates = np.log(np.array(new_rates))
print new_log_rates
[-7.5 -6.5]
[-7.5   -5.875]
[-7.5  -5.25]
[-7.5   -4.625]
/Users/rwest/anaconda/envs/rmg6/lib/python2.7/site-packages/ipykernel/__main__.py:8: RuntimeWarning: divide by zero encountered in log10
[-7.5 -4. ]
[-7.5   -3.375]
[-7.5  -2.75]
[-7.5   -2.125]
[-7.5 -1.5]
[-6.8125 -6.5   ]
[-6.8125 -5.875 ]
[-6.8125 -5.25  ]
[-6.8125 -4.625 ]
[-6.8125 -4.    ]
[-6.8125 -3.375 ]
[-6.8125 -2.75  ]
[-6.8125 -2.125 ]
[-6.8125 -1.5   ]
[-6.125 -6.5  ]
[-6.125 -5.875]
[-6.125 -5.25 ]
[-6.125 -4.625]
[-6.125 -4.   ]
[-6.125 -3.375]
[-6.125 -2.75 ]
[-6.125 -2.125]
[-6.125 -1.5  ]
[-5.4375 -6.5   ]
[-5.4375 -5.875 ]
[-5.4375 -5.25  ]
[-5.4375 -4.625 ]
[-5.4375 -4.    ]
[-5.4375 -3.375 ]
[-5.4375 -2.75  ]
[-5.4375 -2.125 ]
[-5.4375 -1.5   ]
[-4.75 -6.5 ]
[-4.75  -5.875]
[-4.75 -5.25]
[-4.75  -4.625]
[-4.75 -4.  ]
[-4.75  -3.375]
[-4.75 -2.75]
[-4.75  -2.125]
[-4.75 -1.5 ]
[-4.0625 -6.5   ]
[-4.0625 -5.875 ]
[-4.0625 -5.25  ]
[-4.0625 -4.625 ]
[-4.0625 -4.    ]
[-4.0625 -3.375 ]
[-4.0625 -2.75  ]
[-4.0625 -2.125 ]
[-4.0625 -1.5   ]
[-3.375 -6.5  ]
[-3.375 -5.875]
[-3.375 -5.25 ]
[-3.375 -4.625]
[-3.375 -4.   ]
[-3.375 -3.375]
[-3.375 -2.75 ]
[-3.375 -2.125]
[-3.375 -1.5  ]
[-2.6875 -6.5   ]
[-2.6875 -5.875 ]
[-2.6875 -5.25  ]
[-2.6875 -4.625 ]
[-2.6875 -4.    ]
[-2.6875 -3.375 ]
[-2.6875 -2.75  ]
[-2.6875 -2.125 ]
[-2.6875 -1.5   ]
[-2.  -6.5]
[-2.    -5.875]
[-2.   -5.25]
[-2.    -4.625]
[-2. -4.]
[-2.    -3.375]
[-2.   -2.75]
[-2.    -2.125]
[-2.  -1.5]
[  0.           0.          16.10066795  13.42425935  11.42044551
  11.39016056  11.38955889  11.38859707  11.38484183   0.          18.36698786
  16.41404958  13.31346802  11.66201007  11.64809758  11.64878867
  11.64728949  11.38686417  20.45959409  18.8196844   15.70657192
  12.13548437  11.00369217  11.00080002  11.00052683  10.99888818
  10.15159937  20.82801933  17.31749297  13.7705665   10.35646314
   9.76484649   9.763328     9.76514589   9.43654767   8.20317553
  19.01347925  15.40509473  11.76751552   8.77187688   8.48293643
   8.48411869   8.45834356   7.47902287   5.92947147  16.99133029
  13.40788746   9.78603411   7.18899645   7.06467594   7.06621966
   6.56555603   5.12511703   0.          15.01342615  11.41252902
   7.84898947   5.58256359   5.53327411   5.45596382   4.40592296   0.           0.
  12.97732847   9.38092392   6.0756855    4.1616871    4.14427538
   3.60825734   0.           0.           0.          10.83885285
   7.42304491   3.81332715   2.60496671   2.47409188   0.           0.           0.
   0.        ]
In [138]:
new_rate_grid = np.reshape(new_log_rates, (grid_size,grid_size))
plt.imshow(new_rate_grid.T, interpolation='none', origin='lower', extent=extent2, aspect='equal')
for metal, coords in abildpedersen_energies.iteritems():
    color = {'Ag':'w','Au':'w','Cu':'w'}.get(metal,'k')
    plt.plot(coords[0], coords[1], 'o'+color)
    plt.text(coords[0], coords[1], metal, color=color)
plt.xlim(carbon_range)
plt.ylim(oxygen_range)
plt.xlabel('$\Delta E^C$ (eV)')
plt.ylabel('$\Delta E^O$ (eV)')
Out[138]:
<matplotlib.text.Text at 0x1356bae10>
In [68]:
def get_reaction_count(experiment):
    d = directory(*experiment)
    f = os.path.join(d,'chemkin','chem_annotated.inp')
    with open(f) as chemkin:
        r = re.compile('\! Reaction index: Chemkin \#(\d+)')
        for line in chemkin:
            m = r.match(line)
            if m:
                count = int(m.group(1))
    return count
    
get_reaction_count(experiment)
Out[68]:
33
In [134]:
i=35
print experiments[i]
print get_reaction_count(experiments[i])
[-5.4375 -1.5   ]
448
In [143]:
reaction_counts = map(get_reaction_count, experiments)
reaction_counts_grid = np.log10(np.reshape(reaction_counts, (grid_size,grid_size)))
plt.imshow(reaction_counts_grid.T, interpolation='none', origin='lower', extent=extent2, aspect='equal')
for metal, coords in abildpedersen_energies.iteritems():
    color = {'Ag':'w','Au':'w','Cu':'w'}.get(metal,'k')
    color='k'
    plt.plot(coords[0], coords[1], 'o'+color)
    plt.text(coords[0], coords[1], metal, color=color)
plt.xlim(carbon_range2)
plt.ylim(oxygen_range2)
plt.xlabel('$\Delta E^C$ (eV)')
plt.ylabel('$\Delta E^O$ (eV)')

for e,n in zip(experiments,reaction_counts):
    plt.text(e[0],e[1],n,color='w',ha='center', va='center')
#plt.colorbar()
In [146]:
reaction_counts_grid = np.reshape(reaction_counts, (grid_size,grid_size))
ax = sns.heatmap(reaction_counts_grid.T[::-1,:], annot=True, fmt='d', square=True)

STOP HERE.

stuff below is left over from old notebook

In [45]:
raise NotImplementedError("Stop here.")
---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
<ipython-input-45-f3b1bc78107e> in <module>()
----> 1 raise NotImplementedError("Stop here.")

NotImplementedError: Stop here.

Model generation: with RMG-Cat reaction families

Now lets look at the version that generates a mechanism by applying reaction families. First, inspect how the input file differs from the one above.

In [ ]:
%%bash
python $RMGpy/rmg.py ch4_co2_families/input.py > /dev/null
tail -n12 ch4_co2_families/RMG.log

Data processing

First, we make a few plots of the new model. Mostly just to show how it can be done, and to see what the results look like. These aren't very pretty as they're not going in the manuscript.

In [ ]:
data2 = get_pandas_data('ch4_co2_families')
rename_columns(data2)
data2.columns
In [ ]:
get_last_csv_file('ch4_co2_families')
In [ ]:
data2[['CH4', 'CO2']].plot.line()
data2[['CO', 'H2']].plot.line()
data2[['H2O']].plot.line()
In [ ]:
print "All species"
data2.plot.area()
plt.show()
In [ ]:
print "Significant species (those that exceed 0.001 mol at some point)"
significant = [n for n in data2.columns if(data2[n]>0.001).any()]
with sns.color_palette("hls", len(significant)):
    data2[significant].plot.area(legend='reverse')
In [ ]:
surface = [n for n in data2.columns if 'X' in n and n!='X' and (data2[n]>1e-6).any()]
print "The {} surface species that exceed 1e-6 mol at some point".format(len(surface))
with sns.color_palette('Set3',len(surface)):
    data2[surface].plot.area(legend='reverse')
In [ ]:
species_names = data2.columns
gas_phase = [n for n in species_names if 'X' not in n and (data2[n]>0).any()]
surface_phase = [n for n in species_names if 'X' in n]
surface_phase.remove('X')
print "Total moles of gas"
data2[gas_phase].sum(axis=1).plot()
plt.show()
print "All gas phase species"
data2[gas_phase].plot.line()
plt.show()
print "All surface species"
data2[surface_phase].plot.line()
plt.show()
In [ ]:
print "A comparison of the two reactants across the two models"
ax = plt.subplot()
data1[['CH4', 'CO2']].plot.line(ax=ax) # seed
data2[['CH4', 'CO2']].plot.line(ax=ax, style=':') # families
plt.show()

Model comparison

Now we will start making some prettier plots for the manuscript, comparing the two models

In [ ]:
plt.rcParams.update({'mathtext.default':  'regular' }) # make the LaTeX subscripts (eg. CH4) use regular font
sns.set_context("paper", font_scale=1, rc={"lines.linewidth": 1.5}) # Tweak the font size and default line widths

def comparison_plot(subplot_axis=None, species='CH4'):
    label = re.sub('(\d)',r'$_\1$',species)
    ax1 = subplot_axis or plt.subplot()
    
    plt.locator_params(nbins=4) # fewer tick marks
    
    ax1.plot(0, data1[species].iloc[0], 'ko')
    ax1.plot(data1.index, data1[species],'k--', linewidth=2.5)
    ax1.plot(data2.index, data2[species],'r-')
    plt.xlim(0,1)
    ax1.set_ylabel('moles')
    ax1.set_xlabel('time (s)')
    
    # Legend
    dummy = matplotlib.patches.Patch(alpha=0)
    plt.legend([dummy],[label], loc='best', fontsize=12)
    
    plt.tight_layout()

fig = plt.figure(figsize=(2.5,2.5))
ax1 = comparison_plot()
In [ ]:
fig = plt.figure(figsize=(7,4))
for n, species in enumerate('CH4 CO H2O CO2 H2'.split()):
    ax = plt.subplot(2,3,n+1)
    comparison_plot(ax, species)

ax = plt.subplot(2,3,6)
ax.axis('off')
red = matplotlib.lines.Line2D([], [], color='r', label='RMG-Cat')
dash = matplotlib.lines.Line2D([], [], color='k', linestyle='--', linewidth=2.5, label='Delgado et al.')

plt.legend(handles=[red, dash], loc='best',fontsize=12)
plt.tight_layout()
plt.savefig('Multi-panel gas comparison.pdf', bbox_inches='tight')
In [ ]:
def extras_plot(subplot_axis=None, species_list=['C2H6','CH2O'], label_positions=None):
    """
    Plot the requested species on one plot,
    with just the 'data2' (RMG-Cat) values.
    Useful for species not in the Delgado model.
    """
    ax1 = subplot_axis or plt.subplot()
    plt.locator_params(nbins=4) # fewer tick marks
    plt.ticklabel_format(style='sci', axis='y', scilimits=(3,3))
    
    for i,species in enumerate(species_list):
        label = re.sub('(\d)',r'$_\1$',species)
        ax1.semilogy(data2.index, data2[species],'-', label=label)
        plt.ylim(ymin=1e-9)
        plt.xlim(0,1)
        ax1.set_ylabel('moles')
        ax1.set_xlabel('time (s)')
        # Manually constructed label
        x = label_positions[species] if label_positions and species in label_positions else 0.5
        y = data2[x:][species].iloc[0]
        ax1.text(x,y,label,
                 verticalalignment='bottom',
                 color=sns.color_palette()[i],
                 fontsize=10)
        plt.tight_layout()

# test it
fig = plt.figure(figsize=(4,4))
with sns.color_palette("Dark2", 3):
    extras_plot(species_list='C2H6 CH2O C3H8'.split())
In [ ]:
# Everything that exceeds some lower threshold
gas_phase = [n for n in species_names if 'X' not in n and (data2[n]>1.1e-9).any()]
# Remove things already in the comparison plots
gas_phase = [n for n in gas_phase if n not in 'CH4 CO H2O CO2 H2 N2'.split()]
print "Extra gas phase species of note are {}.".format(", ".join(gas_phase))
fig = plt.figure(figsize=(3.25,3.25))
with sns.color_palette("Dark2", len(gas_phase)):
    extras_plot(species_list=gas_phase, 
               label_positions={'C2H5':0.15, 'C3H8':0.1, 'CH3OH':0.7, 'H':0.3})
plt.savefig('Gas extra species (many).pdf',bbox_inches='tight')
In [ ]:
def multi_comparison_plot(subplot_axis=None, species_list=['X', 'HX'], label_positions=None):
    """
    Plot many surface species AND their comparison (if there is one) from the Delgado model
    on a semilog plot. 
    """
    ax1 = subplot_axis or plt.subplot()
    plt.locator_params(nbins=4) # fewer tick marks
    for i,species in enumerate(species_list):
        assert 'X' in species, "For surface species only (y axis is normalized for fractional coverage)."
        label = re.sub('(\d)',r'$_\1$',species)
        label = label.replace('X','*')
        if label=='*': label = 'vacant'
        if species in data1.columns:
            sites = data1['X'].iloc[0]
            ax1.semilogy(data1.index, data1[species]/sites,'k--', linewidth=2.5)
        sites = data2['X'].iloc[0]
        ax1.semilogy(data2.index, data2[species]/sites,'-')
        plt.xlim(0,1)
        plt.ylim(ymin=1e-6)
        ax1.set_ylabel('site fraction')
        ax1.set_xlabel('time (s)')
        ## Manually constructed label
        x = label_positions[species] if label_positions and species in label_positions else 0.5
        y = data2[x:][species].iloc[0] / sites
        ax1.text(x,y,label, 
                 fontsize=10,
                 verticalalignment='bottom',
                 color=sns.color_palette()[i])

# test it
multi_comparison_plot()
In [ ]:
species_names = data2.columns
gas_phase = [n for n in species_names if 'X' not in n and (data2[n]>0).any()]

sites = data2['X'].iloc[0] # initial number of suface sites
# get only adsorbates that exceed a site fraction of 1e-6
surface_phase = [n for n in species_names if 'X' in n and (data2[n]>sites*1e-6).any()]
surface_phase
In [ ]:
label_positions = {'HX':0.4, 'COX':0.6, 'CHX':0.7, 'OX':0.4, 
                   'CHOX':0.6, 'CX':0.3, 'CH2X':0.05,
                   'C2HOX':0.8, 'CHO2X':0.7}
fig = plt.figure(figsize=(5,5))
with sns.color_palette("Dark2", len(surface_phase)):
    ax1 = multi_comparison_plot(species_list=surface_phase, label_positions=label_positions)
plt.tight_layout()
plt.savefig('Surface comparison semilog.pdf', bbox_inches='tight')
In [ ]:
def multi_comparison_loglogplot(subplot_axis=None, species_list=['X', 'HX'], label_positions=None):
    """
    Plot many things AND their comparison (if there is one) from the Delgado model
    on a log-log plot
    """
    ax1 = subplot_axis or plt.subplot()
    for i,species in enumerate(species_list):
        label = re.sub('(\d)',r'$_\1$',species)
        label = label.replace('X','*')
        if label=='*': label='vacant'
        if species in data1.columns:
            sites = data1['X'].iloc[0]
            ax1.loglog(data1.index, data1[species]/sites,'k--', linewidth=2.5)
        sites = data2['X'].iloc[0]
        ax1.loglog(data2.index, data2[species]/sites,'-')
        plt.xlim(1e-6,1)
        plt.ylim(ymin=1e-6)
        ax1.set_ylabel('site fraction')
        ax1.set_xlabel('time (s)')
        ## Manually constructed label
        x = label_positions[species] if label_positions and species in label_positions else 3.0
        x = 10**(-1*x)
        y = data2[x:][species].iloc[0] / sites
        if x==1: x=1.3 # move them right just a little
        ax1.text(x,y,label, fontsize=10,
                 verticalalignment='center' if x==1.3 else 'bottom',
                 color=sns.color_palette()[i])  

label_positions = {'HX':0.5, 'COX':0, 'CX':2, 'C2HOX':0,
                  'CH3CXO':0, 'CH3CX':0, 'OX':5, 'CHO2X':0,
                  'CH2X':2.5}
fig = plt.figure(figsize=(5,4))
with sns.color_palette("Dark2", len(surface_phase)):
    ax1 = multi_comparison_loglogplot(species_list=surface_phase, label_positions=label_positions)
plt.tight_layout()
plt.savefig('Surface comparison loglog.pdf',bbox_inches='tight')

Effect of tolerance

Now we investigate the effect of gradually decreasing (tightening) the tolerance

In [ ]:
# First, make a series of input files in separate directories

base_directory = 'ch4_co2_tolerances'

with open(os.path.join(base_directory, 'input.template.py')) as infile:
    input_file = infile.read()

def directory(i):
    return os.path.join(base_directory, "ch4_co2_tolerance_m{}".format(i))

for i in range(9):
    tolerance = 0.1**i
    tolerance_string = 'toleranceMoveToCore={:.1e},'.format(tolerance)
    print tolerance_string
    input_file = re.sub('toleranceMoveToCore\s*=\s*(.*?),', tolerance_string, input_file)
    os.path.exists(directory(i)) or  os.makedirs(directory(i))
    with open(os.path.join(directory(i), 'input.py'), 'w') as outfile:
        outfile.write(input_file)
    print "Saved to {}/input.py".format(directory(i))
In [ ]:
# Now run all the jobs
# Don't execute this cell unless you have a while to wait.
import subprocess
import sys
for i in range(9):
    print "Attempting to run job {} in directory {}".format(i, directory(i))
    try:
        retcode = subprocess.call("python $RMGpy/rmg.py {}/input.py".format(directory(i)), shell=True)
        if retcode < 0:
            print >>sys.stderr, "Process was terminated by signal", -retcode
        elif retcode > 0:
            print >>sys.stderr, "Process returned", retcode
        else:
            print "Success"
    except OSError as e:
        print >>sys.stderr, "Execution failed:", e
In [ ]:
# Now read the ends of the log files and extract the mechanism sizes.
epsilon = []
core_species = []
core_rxns = []
edge_species = []
edge_rxns = []

for i in range(9):
    dirname = directory(i)
    print "reading from ", directory(i)
    last_lines = subprocess.check_output(['tail', '-n','6', os.path.join(directory(i),'RMG.log')])

    match = re.search('The final model core has (\d+) species and (\d+) reactions', last_lines)
    if match is None: 
        print "Trouble with {}/RMG.log:\n{}".format(directory(i),last_lines)
    core_species.append(int(match.group(1)))
    core_rxns.append(int(match.group(2)))
    match = re.search('The final model edge has (\d+) species and (\d+) reactions', last_lines)
    if match is None: 
        print "Trouble with {}/RMG.log:\n{}".format(directory(i),last_lines)
    edge_species.append(int(match.group(1)))
    edge_rxns.append(int(match.group(2)))
    epsilon.append( 0.1**i )

#remove the four inerts, N2, Ar, Ne, and He that RMG automatically adds
core_species = np.array(core_species) - 4
edge_species = np.array(edge_species) - 4
In [ ]:
# Now count the reactions by type
Deutschmann = []
abstraction = []
adsorption = []
dissociation = []
for i in range(9):
    ckfilepath = os.path.join(directory(i), 'chemkin', 'chem_annotated.inp')
    re_library = re.compile('! Library reaction: (.*)')
    re_family = re.compile('! Template reaction: (.*)')
    from collections import Counter
    counts = Counter()
    with open(ckfilepath) as ckfile:
        for line in ckfile:
            match = re_family.match(line) or re_library.match(line)
            if match is None:
                continue
            source = match.group(1)
            counts[source] += 1
    counts

    Deutschmann.append(counts['Deutschmann_Ni'])
    abstraction.append(counts['Surface_Abstraction'])
    adsorption.append(counts['Surface_Adsorption_Dissociative'] + counts['Surface_Adsorption_Single'])
    dissociation.append(counts['Surface_Dissociation'])
print 'Deutschmann',Deutschmann
print 'abstraction',abstraction
print 'adsorption',adsorption
print 'dissociation',dissociation

Deutschmann = np.array(Deutschmann)
abstraction = np.array(abstraction)
adsorption = np.array(adsorption)
dissociation = np.array(dissociation)
total = Deutschmann + abstraction + adsorption + dissociation
assert (total == np.array(core_rxns)).all(), "Sum of counters doesn't equal core_rxns from above"
In [ ]:
# Now plot the figure
fig = plt.figure()
fig.set_size_inches(7,3.5) #set size, in inches
gs = matplotlib.gridspec.GridSpec(1, 3)
ax0 = plt.subplot(gs[0])
ax1 = plt.subplot(gs[1])
ax2 = plt.subplot(gs[2])

ax0.loglog(epsilon, edge_species, 'r^-', label='edge')
ax0.loglog(epsilon, core_species, 'bo-', label='core')

ax1.loglog(epsilon, edge_rxns, 'r^-', label='edge')
ax1.loglog(epsilon, core_rxns, 'bo-', label='core')

#ax2.semilogx(epsilon, Deutschmann, 'k:', label='Deutschmann')
#ax2.semilogx(epsilon, abstraction, 'b-', label='abstraction')
#ax2.semilogx(epsilon, dissociation, 'g--', label='dissociation')
#ax2.semilogx(epsilon, adsorption, 'r-.', label='adsorption')

stacks = ax2.stackplot(epsilon, Deutschmann, adsorption, dissociation, abstraction,
             labels='Deutcschmann Adsorption Dissociation  Abstraction'.split(),
             colors=sns.color_palette("Set2",4)
             )
hatches = ['', '///', '|||', '---']
for i, patch in enumerate(stacks):
    patch.set_hatch(hatches[i])
ax2.set_xscale('log')


#ax0.set_ylim([20,1000])
#ax1.set_ylim([30,2000])
#ax2.set_ylim([0,300])

ax0.set_xlabel('error tolerance, $\epsilon$', fontsize=11)
ax1.set_xlabel('error tolerance, $\epsilon$', fontsize=11)
ax2.set_xlabel('error tolerance, $\epsilon$', fontsize=11)

#ax0.set_ylabel('species')
#ax1.set_ylabel('reactions')
#ax2.set_ylabel('reactions in core')

ax0.set_title('no. of species', fontsize=11)
ax1.set_title('no. of reactions', fontsize=11)
ax2.set_title('no. of core reactions', fontsize=11)

ax0.legend(loc='upper left', bbox_to_anchor=[0.0, 0.9], fontsize=10, frameon=False)
ax1.legend(loc='upper left', bbox_to_anchor=[0.0, 0.9], fontsize=10, frameon=False)
# For the stacked plot we want to plot the legend in reverse order 
# so the bottom patch is on the bottom of the legend
handles, labels = ax2.get_legend_handles_labels()
ax2.legend(handles[::-1], labels[::-1],
          loc='upper left', bbox_to_anchor=[0.0, 0.9], fontsize=10, frameon=False)


ax0.set_xticks([1.0E-8, 1.0E-4,  1.0])
ax1.set_xticks([1.0E-8, 1.0E-4,  1.0])
ax2.set_xticks([1.0E-8, 1.0E-4,  1.0])
ax0.tick_params(labelsize=10)
ax1.tick_params(labelsize=10)
ax2.tick_params(labelsize=10)

ax0.invert_xaxis()
ax1.invert_xaxis()
ax2.invert_xaxis()

ax0.text(0.1, 0.9, '(a)', fontsize=10,
verticalalignment='bottom', horizontalalignment='left',
transform=ax0.transAxes)
ax1.text(0.1, 0.9, '(b)', fontsize=10,
verticalalignment='bottom', horizontalalignment='left',
transform=ax1.transAxes)
ax2.text(0.1, 0.9, '(c)', fontsize=10,
verticalalignment='bottom', horizontalalignment='left',
transform=ax2.transAxes)

fig.tight_layout()
fig.savefig('mechanism_size.pdf', bbox_inches='tight')
In [ ]: